AM, FM SNR calculation in Matlab

amsnr.m

%amsnr: SNR in AM systems
%
%This code contains functions to generate a message signal, amplitude
%modulate a carrier and detect the message using envelope detection. 
%
%Homework: Write code to perform SNR analysis as discussed in homework
%5. Additional information is available in the comments below.
%
% Last modified: 2009-April-07
%
% (c) P. Ramuhalli. All Rights Reserved.
%

%Note the functional approach to solving this problem: We create and use
%sub-routines for each of the major tasks and call these subroutines in the
%main function.

function amsnr
close all, clear all
%This is the main function. This will call the modulation and demodulation
%routines and perform the SNR analysis here.
%Note that we do not have any inputs or outputs for this function.
mu = [0.1 0.3 0.5 0.8 1];  %Modulation index
%Close all open figures so we can start with a clean slate.

for outer = 1:length(mu)
%STEP 1:
%Create a message (sinusoid)
[msg,tim,B,Fs] = generate_msg;%See function definition for generate_msg below for an explanation of the variables.
N=numel(msg);
%STEP 2:
%Now, create a modulated signal with no noise. To do so, we first set up
%the required variables.

A = sqrt(2); %Carrier amplitude for unity power carrier
fc = 1e2; %Carrier frequency in Hz
f = 0:Fs/2;
%Now call the AM modulation function
[sam] = AMmod(msg,tim,A,fc,mu(outer)); %See function definition for AMmod below for an explanation of the variables.
%In the line above, sam is the AM signal.

%Now, add noise with given noise power. It may be better to put this
%section in a loop to cover all desired SNR values. We want to compute
%the received SNR, predetection SNR and post-detection SNR.

%Create a list of desired SNR values
SNRvals = [1,5,10,15,20,25,30,35];
Nfact=10.^(SNRvals./10); %factor by which powers will differ

%Get power of original modulated signal
%Si=Power(sam,N,f,Fs);
Si=(A^2+mean((mu(outer)*msg).^2))/2;
mdet_clean=env_demod(sam,B,Fs);
mhat_clean=DCBlock(mdet_clean);
Namp = sqrt(Si./Nfact); %finds new power of noise and converts to ampl.
%In a loop, calculate the different SNR values
for jj = 1:length(SNRvals)

    %For each desired SNR value, DO: (YOU NEED TO ADD CODE TO DO THIS!)
    
    %STEP 3: Calculate AWGN noise power for AM and add noise to signal. 
    %Note that the noise power is relative to signal power, and SNR is in dB
    noise(jj,:)=Namp(jj).*randn(size(sam));
    Sair(jj,:)=sam+noise(jj,:);
    %Sair=awgn(sam,SNRvals(jj),'measured');
    
    %STEP 4: Compute SNR at receiver front-end. 
    %This will be close to the desired
    %SNR, but not the same due to numerical errors and small number of
    %random samples.
    Ni(jj)=mean(noise(jj,:).^2);
    SNRair(jj)=10*log10(Si/Ni(jj));

    %STEP 5: Calculate Predetection SNR
    %First, send the noisy signal throught the BPF
    %We also transmit only the signal through the BPF
    %Now calculate SNR. Noise power is calculated using the difference between the signal and noisy signal.
    Spre(jj,:) = BPF(Sair(jj,:),N,Fs,B,fc); 
    Sprepwr(jj) = Power(Spre(jj,:),N,f,Fs);
    Nprepwr(jj) = mean((Spre(jj,:)-Sair(jj,:)).^2);
    SNRpre(jj) = 10*log10(Sprepwr(jj)./Nprepwr(jj));
    
    
    %STEP 6:Envelope Detection and Postdetection SNR
    %Again, we demodulate the noisy BPF filtered signal first.
    %Now, we demodulate the clean BPF filtered signal.
    %Finally, we DC Block both the detected signals above
    %Now calculate postdetection SNR.
    mdet_air(jj,:)=env_demod(Sair(jj,:),B,Fs);
    mhat_air(jj,:)=DCBlock(mdet_air(jj,:));
    So(jj) = mean(mhat_air(jj,:).^2);
    SNRpost(jj) = 10*log10(Si./So(jj));
    
    
    %We repeat this process in a loop to cover all desired SNR values. The
    %calculated SNR values are held in an array for plotting.
end

%Plot pre-detection vs post detection SNR values (ADD CODE TO DO THIS!)


subplot(5,1,outer)
stem(SNRpre,SNRpost)
title({['SNR for mu=',num2str(mu(outer))]})
xlabel('SNRpre'),ylabel('SNRpost')
end

function SigPwr = Power(sig,N,f,Fs)
   
SigPwr = sum(pwelch(sig,ones(N,1),0,N,Fs,'twosided')./Fs);
%SigPwr = sum(psd(sig));
%SigPwr=sum(abs(fft(sig)).^2/N);
end

function  [m,t,B,Fs] = generate_msg
%This subroutine generates a sinusoidal message. 
%Inputs: None
%Outputs:
%   m: Message
%   t: time axis
%   B: BW of message
%   Fs: Sampling frequency
%

%Set basic parameters
fm=2; %Carrier frequency in Hz
Fs = 1e3; %Sampling frequency in Hz
B = 10; %Approximate BW of the message
%changed B to 10Hz from 

%Create message signal
t = 0:1/Fs:2; %Time axis
m = sin(2*pi*(fm.*t)); %message signal

end

function  [sam] = AMmod(m,t,A,fc,mu)
%This subroutine generates an AM modulated signal. We assume that the
%modulating signal is a square wave.
%Inputs:
%   m: Message
%   t: time axis
%`  A: Carrier amplitude
%   fc: Carrier frequency
%   mu: Modulation index
%Outputs:
%   sam: AM modulated signal
%

sam = A*(1+(mu*m)).*cos((2*pi*fc).*t); %AM modulated signal
end


function [m]=env_demod(s,B,Fs)
%Envelope demodulation is equivalent to rectification followed by low pass
%filtering
%Inputs:
%   s: Amplitude modulated signal (could also be slope-detected FM signal)
%   B: BW of message
%   Fs: Sampling frequency
%Outputs:
%   m: Demodulated message
%

%Rectify the AM signal
s1 = s.*(s>0); %Half wave rectification

%Now design a low pass filter with passband equal to B
%We first compute the normalized frequency between 0 and 1, with 1
%corresponding to half the sample rate, Fs/2. We assume that B < Fs/2, but
%do not check for this!

N    = 256;       % Filter Order
wc = (B/2)/(Fs/2); %Cutoff frequency
win = hamming(N+1); %Use Hamming window
h = fir1(N,wc, 'low', win, 'scale'); %Design FIR filter
%Now filter:
s1 = [s1(N:-1:1),s1];
m = filter(h,1,s1);
m = m(N+1:length(m));
%We will have some initial distortion due to the convolution operation
%being over a finite time period. But we ignore this in our presentation.

end


function BPFd = BPF(s1,N,Fs,B,fc)
N    = 256;       % Filter Order
wc = [fc-B fc+B]./(Fs/2); %Cutoff frequency
win = hamming(N+1); %Use Hamming window
h = fir1(N,wc, 'low', win, 'scale'); %Design FIR filter
%Now filter:
s1 = [s1(N:-1:1),s1];
BPFd = filter(h,1,s1);
BPFd = BPFd(N+1:length(BPFd));
%We will have some initial distortion due to the convolution operation
%being over a finite time period. But we ignore this in our presentation.

end


function mhat = DCBlock(m)
% DC blocking is equivalent to mean removal
%Inputs:
%   m: Demodulated message
%Outputs:
%   mhat: DC shifted message
%

%Subtract time-mean from signal.
mhat = m-mean(m);
end

fmsnr.m

%fmsnr: SNR in FM systems
%
%This code contains functions to generate a message signal, amplitude
%modulate a carrier and detect the message using envelope detection. 
%
%Homework: Write code to perform SNR analysis as discussed in homework
%5. Additional information is available in the comments below.
%
% Last modified: 2009-April-07
%
% (c) P. Ramuhalli. All Rights Reserved.
%

%Note the functional approach to solving this problem: We create and use
%sub-routines for each of the major tasks and call these subroutines in the
%main function.

function fmsnr;
%This is the main function. This will call the modulation and demodulation
%routines and perform the SNR analysis here.
%Note that we do not have any inputs or outputs for this function.

%Close all open figures so we can start with a clean slate.
close all, clear all


%STEP 1:
%Create a message (sinusoid)
[msg,tim,Fs] = generate_msg;%See function definition for generate_msg below for an explanation of the variables.
N=numel(msg);
f = 0:Fs/2;
%STEP 2:
%Now, create a modulated signal with no noise. To do so, we first set up
%the required variables.
kf = 2*pi; %Modulation index
A = sqrt(2); %Carrier amplitude for unity power carrier
fc = 1e3; %Carrier frequency in Hz
%Now call the AM modulation function
[sem,B] = FMmod(msg,tim,A,fc,kf); %See function definition for FMmod below for an explanation of the variables.
%In the line above, sem is the FM signal.

%Check signal power
h_sig = spectrum.welch;                  % Create a Welch spectral estimator.
PSD_sem = psd(h_sig,sem,'Fs',Fs);        %Compute the spectrum
pwr_sem = avgpower(PSD_sem);             %Use the avgpower () function from MATLAB to compute the power

%Now, add noise with given noise power. It may be better to put this
%section in a loop to cover all desired SNR values. We want to compute
%the received SNR, predetection SNR and post-detection SNR.

%Create a list of desired SNR values
SNRvals = [1,5,10,15,20,25,30, 35];
Nfact=10.^(SNRvals./10); %factor by which powers will differ
Si=A^2/2;
Namp = sqrt(Si./Nfact);
%In a loop, calculate the different SNR values
for jj = 1:length(SNRvals)

    %For each desired SNR value, DO: (YOU NEED TO ADD CODE TO DO THIS!)
    
    %STEP 3: Calculate AWGN noise power for FM. 
    %Note that the noise power is relative to signal power, and SNR is in dB
     noise(jj,:)=Namp(jj).*randn(size(sem));
    Sair(jj,:)=sem+noise(jj,:);
    
    %STEP 4: Compute SNR at receiver front-end. 
    %This will be close to the desired
    %SNR, but not the same due to numerical errors and small number of
    %random samples.
        Ni(jj)=mean(noise(jj,:).^2);
    SNRair(jj)=10*log10(Si/Ni(jj));

    %STEP 5: Calculate Predetection SNR
    %First, send the noisy signal throught the BPF
    %We also transmit only the signal through the BPF
    %Now calculate SNR. Noise power is calculated using the difference between the signal and noisy signal.
        Spre(jj,:) = BPF(Sair(jj,:),N,Fs,B,fc); 
    Sprepwr(jj) = Power(Spre(jj,:),N,f,Fs);
    Nprepwr(jj) = mean((Spre(jj,:)-Sair(jj,:)).^2);
    SNRpre(jj) = 10*log10(Sprepwr(jj)./Nprepwr(jj));
    
    %STEP 6:Envelope Detection and Postdetection SNR
    %Again, we demodulate the noisy BPF filtered signal first.
    %The first step is time derivative. This is followed by
    %envelope-detection and DC blocking.
    mdet_air(jj,:)=env_demod(Sair(jj,:),B,Fs);
    mhat_air(jj,:)=DCBlock(mdet_air(jj,:));
    So(jj) = mean(mhat_air(jj,:).^2);
    SNRpost(jj) = 10*log10(Si./So(jj));
        
    %Now, we demodulate the clean BPF filtered signal.
    %Now calculate postdetection SNR.
    
    %We repeat this process in a loop to cover all desired SNR values. The
    %calculated SNR values are held in an array for plotting.
    stem(SNRpre,SNRpost)
title('SNR FM')
xlabel('SNRpre'),ylabel('SNRpost')
end

%Plot pre-detection vs post detection SNR values (ADD CODE TO DO THIS!)


function  [m,t,Fs] = generate_msg;
%This subroutine generates a sinusoidal message. 
%Inputs: None
%Outputs:
%   m: Message
%   t: time axis
%   Fs: Sampling frequency
%

%Set basic parameters
fm=1; %Carrier frequency in Hz
Fs = 5e4; %Sampling frequency in Hz

%Create message signal
t = 0:1/Fs:2; %Time axis
m = sin(2*pi*(fm.*t)); %message signal

end

function  [sem,B] = FMmod(m,t,A,fc,kf);
%This subroutine generates an FM modulated signal. 
%Inputs: 
%   m: Message
%   t: time axis
%`  A: Carrier amplitude
%   fc: Carrier frequency
%   kf: Modulation constant
%Outputs: 
%   sem: FM modulated signal

%Generate FM
f_fm=fc+(kf*m); %Instantaneous frequency
sem = A*cos(2*pi*(f_fm.*t)); %FM modulated signal

B = 200;
end


function dsem = TimeDeriv(sem);
%This function performs the differentiation operation
%Inputs: 
%   sem: FM signal
%Outputs: 
%   dsem: Derivative of FM signal (slope detected FM signal)
%

%Compute Derivative
dsem = [0,diff(sem)];
end

function [m]=env_demod(s,B,Fs);
%Envelope demodulation is equivalent to rectification followed by low pass
%filtering
%Inputs:
%   s: Amplitude modulated signal (could also be slope-detected FM signal)
%   B: BW of message
%   Fs: Sampling frequency
%Outputs:
%   m: Demodulated message
%

%Rectify the AM signal
s1 = s.*(s>0); %Half wave rectification

%Now design a low pass filter with passband equal to B
%We first compute the normalized frequency between 0 and 1, with 1
%corresponding to half the sample rate, Fs/2. We assume that B < Fs/2, but
%do not check for this!

N    = 256;       % Filter Order
wc = (B/2)/(Fs/2); %Cutoff frequency
win = hamming(N+1); %Use Hamming window
h = fir1(N,wc, 'low', win, 'scale'); %Design FIR filter
%Now filter:
s1 = [s1(N:-1:1),s1];
m = filter(h,1,s1);
m = m(N+1:length(m));
%We will have some initial distortion due to the convolution operation
%being over a finite time period. But we ignore this in our presentation.

end

function mhat = DCBlock(m);
% DC blocking is equivalent to mean removal
%Inputs:
%   m: Demodulated message
%Outputs:
%   mhat: DC shifted message
%

%Subtract time-mean from signal.
mhat = m-mean(m);
end

function SigPwr = Power(sig,N,f,Fs)
   
SigPwr = sum(pwelch(sig,ones(N,1),0,N,Fs,'twosided')./Fs);
%SigPwr = sum(psd(sig));
%SigPwr=sum(abs(fft(sig)).^2/N);
end


function BPFd = BPF(s1,N,Fs,B,fc)
N    = 256;       % Filter Order
wc = [fc-B fc+B]./(Fs/2); %Cutoff frequency
win = hamming(N+1); %Use Hamming window
h = fir1(N,wc, 'low', win, 'scale'); %Design FIR filter
%Now filter:
s1 = [s1(N:-1:1),s1];
BPFd = filter(h,1,s1);
BPFd = BPFd(N+1:length(BPFd));
%We will have some initial distortion due to the convolution operation
%being over a finite time period. But we ignore this in our presentation.

end