Search code examples
matlabplotmean-square-error

MATLAB: : Mean square error vs SNR plot


I am having difficulty in understanding the logic behind generating a plot of SNR (db) vs MSE. Different Signal to Noise Ratio (SNR) is created by varying the noise power . The formula of MSE is averaged over T independent runs.

For each SNR, I generate NEval = 10 time series. How do I correctly plot a graph of SNR vs MSE when SNR is in the range = [0:5:50]? Below is the pseudo code.

 N = 100; %Number_data_points
NEval = 10; %Number_of_different_Signals
Snr = [0:5:50];
T = 1000; %Number of independent runs

MSE = [1];
for I = 1:T
 for snr = 1: length(Snr)
   for expt = 1:NEval
    %generate signal
      w0=0.001;  phi=rand(1);
    signal = sin(2*pi*[1:N]*w0+phi); 
    % add zero mean Gaussian noise
    noisy_signal = awgn(signal,Snr(snr),'measured');
    % Call Estimation algorithm
    %Calculate error
    end
  end
end

plot(Snr,MSE); %Where and how do I calculate this MSE

Solution

  • As explained here (http://www.mathworks.nl/help/vision/ref/psnr.html) or other similar sources, MSE is simply the mean squared error between the original and corrupted signals. In your notations,

    w0=0.001;
    signal = sin(2*pi*[1:N]*w0);
    MSE = zeros(T*Neval,length(Snr));
    for snr = 1:length(Snr)
        for I = 1:T*Neval     %%note, T and Neval play the same role in your pseudo code
            noisy_signal = awgn(sin(2*pi*[1:N]*w0+rand(1)),Snr(snr),'measured');
            MSE(I,snr) = mean((noisy_signal - signal).^2);
        end
    end
    semilogy(Snr, mean(MSE))  %%to express MSE in the log (dB-like) units
    

    For the case of signals of different lengths:

    w0=0.001;
    Npoints = [250,500,1000];    
    MSE = zeros(T,length(Npoints),length(Snr));
    for snr = 1:length(Snr)
      for ip = 1:length(Npoints) 
        signal = sin(2*pi*[1:Npoints(ip)]*w0);
        for I = 1:T
            noisy_signal = awgn(sin(2*pi*[1:Npoints(ip)]*w0+rand(1)),Snr(snr),'measured');
            MSE(I,ip,snr) = mean((noisy_signal - signal).^2);
        end
      end
    end
    semilogy(Snr, squeeze(mean(mean(MSE,1),2)) )