%%% Illustrate FM modulation and demodulation using triangl function 
%%% Defining message signal as a triangle signal, in real life however, this
%%% could be anything like voice signal etc. that we would want to transmit
function ExampleFM
clear all;
clf;

ts=1e-4; % sampling interval
t=-0.04:ts:0.04;
Ta=0.01;

% Use triangl function to generate message signal
m_sig=triangl((t+0.01)/Ta)-triangl((t-0.01)/Ta);

%Lm_sig=length(m_sig1);
Lfft=length(t);   % defining DFT (or FFT) size
Lfft=2^ceil(log2(Lfft))  % making Lfft a power of 2 since this makes the fft algorithm work fast
M_fre=fftshift(fft(m_sig,Lfft))  % i.e. calculating the frequency domain message signal,
                                 % fft algorithm calculates points from ) to Lfft-1, hence we use fftshift on this
                                 % result to order samples from -Lfft/2 to Lfft/2 -1 
freqm=(-Lfft/2:Lfft/2-1)/(Lfft*ts) % Defining the frequency axis for the frequency domain DSB modulated signal
B_m=100; % bandwidth of the signal is B_m Hz
h=fir1(80,[B_m*ts]) % defining a simple lowpass filter with bandwidth B_m Hz

kf=160*pi;
m_intg=kf*ts*cumsum(m_sig);
s_fm=cos(2*pi*300*t+m_intg)
s_pm=cos(2*pi*300*t+m_sig)
Lfft=length(t);   % defining DFT (or FFT) size
Lfft=2^ceil(log2(Lfft)+1)  % increasing Lfft by factor of 2 

S_fm=fftshift(fft(s_fm,Lfft)); % obtaining frequency domain modulated signal
S_pm=fftshift(fft(s_pm,Lfft)); % obtaining frequency domain modulated signal
freqs=(-Lfft/2:Lfft/2-1)/(Lfft*ts); % Defining the frequency axis for the frequency domain DSB modulated signal


%%% Demodulation 
% Using an ideal low pass filter with bandwidth 200 Hz
s_fmdem=diff([s_fm(1) s_fm])/ts/kf
s_fmrec=s_fmdem.*(s_fmdem>0);
s_dec=filter(h,1,s_fmrec)

  
Trange=[-0.04 0.04 -1.2 1.2] % axis ranges for signal, this specifies the range of axis for the plot, the first two parameters are range limits for x-axis, and last two parameters are for y-axis

Frange=[-600 600 0 300] % axis range for frequency domain plots

figure(1)
subplot(211); m1=plot(t,m_sig)
axis(Trange) % set x-axis and y-axis limits 
xlabel('{\it t}(sec)'); ylabel('{\it m}({\it t})')
title('Message signal')

subplot(212); m2=plot(t,s_dec)
xlabel('{\it t}(sec)'); ylabel('{\it m}_d({\it t})')
title('Demodulated FM signal')

figure(2)
subplot(211); td1=plot(t,s_fm)
axis(Trange) % set x-axis and y-axis limits 
title('FM signal')
xlabel('{\it t}(sec)'); ylabel('{\it s}_{\rm FM}({\it t})')


subplot(212); td1=plot(t,s_pm)
axis(Trange) % set x-axis and y-axis limits 
title('PM signal')
xlabel('{\it t}(sec)'); ylabel('{\it s}_{\rm PM}({\it t})')




figure(3)
subplot(211); fp1=plot(t,s_fmdem)
%axis(Trange) % set x-axis and y-axis limits 
xlabel('{\it t}(sec)'); ylabel('{\it d s}_{\rm FM}({\it t})')
title('FM Derivative')

subplot(212); fp2=plot(t,s_fmrec)
%axis(Trange) % set x-axis and y-axis limits 
xlabel('{\it t}(sec)'); %ylabel('{\it d}_{\rm PM}({\it t})')
title('rectified FM Derivative')



figure(4)

subplot(211); fd1=plot(freqs,abs(S_fm))
axis(Frange) % set x-axis and y-axis limits 
xlabel('{\it f}(Hz)'); ylabel('{\it S}_{\rm FM}({\it f})')
title('FM Amplitude Spectrum')

subplot(212); fd2=plot(freqs,abs(S_pm))
axis(Frange) % set x-axis and y-axis limits 
xlabel('{\it f}(Hz)'); ylabel('{\it S}_{\rm PM}({\it f})')
title('PM Amplitude Spectrum')





%%% Defining triangl function  used in above code
% triangl(t)=1-|t| , if |t|<1
% triangl(t)=0 , if |t|>1


function y = triangl(t)
y=(1-abs(t)).*(t>=-1).*(t<1); % i.e. setting y to 1 -|t|  if  |t|<1 and to 0 if not
%end

%%% example usage
% t=-5:.1:5
% y=triangl(t)
% stem(t,y)