Wednesday, January 8, 2014

Convergence to signal through repeated sampling to average out white-noise in signals (MATLAB)

Normally a filter such as Guassian filter requires that the user knows that frequency they want to look over. However, if a signal is given to a user without telling him/her where to look for the desired signal frequency, the Guassian filter approach may not work, in this case, we can make use of the idea that white noises are always out of phases and random. In this case, we can repeatedly sample signals containing white noise and average them in the frequency domain, potentially cancelling out the white noise

The MATLAB code below creates a noisy signal (un in time domain and utn in frequency domain) from the sech(t) function (u in time domain and ut in frequency domain

clear all; close all; clc;

T=30;
n=512;

t2=linspace(-T/2, T/2, n+1);

t=t2(1:n);

k=(2*pi/T)*[0:n/2-1 -n/2:-1];
ks=fftshift(k);

u=sech(t);

ut=fft(u);
noise=30;
utn=ut+noise*(randn(1, n)+i*randn(1, n));
un=ifft(utn);

subplot(2, 1, 1); plot(t, u, 'k', t, un);
subplot(2, 1, 2); plot(ks, abs(fftshift(ut)), 'k', ks, abs(fftshift(utn)));

Below shows the noisy signal and the original signal in both time and frequency domain

Image and video hosting by TinyPic

The following MATLAB code creates 30 sample noisy signals then take the average of the signals in the frequency domain, which is then converted back to the time domain using ifft

clear all; close all; clc;

T=30;
n=512;

t2=linspace(-T/2, T/2, n+1);

t=t2(1:n);

k=(2*pi/T)*[0:n/2-1 -n/2:-1];
ks=fftshift(k);

u=sech(t);

ut=fft(u);
noise=30;
avg=zeros(1, n);
realizations=30;
for j=1:realizations % 30 sample signals in time and frequency domain, each with periodic domain of T=30
    utn=ut+noise*(randn(1, n)+i*randn(1, n));
    un=ifft(utn);
    avg=avg+utn;
end

unf=ifft(avg);
avg=abs(fftshift(avg))/realizations;

subplot(2, 1, 1); plot(t, u, 'r', t, un, 'c', t, unf, 'k');
subplot(2, 1, 2); plot(ks, abs(fftshift(ut)), 'r', ks, abs(fftshift(utn)), 'c', ks, avg, 'k');

The figure below shows the filtered results using the convergence via averaging of signal samples in frequency domain

Image and video hosting by TinyPic

The MATLAB code below further shows an animation of convergence as the number of samples taken is increased.

clear all; close all; clc;

T=30;
n=512;

t2=linspace(-T/2, T/2, n+1);

t=t2(1:n);

k=(2*pi/T)*[0:n/2-1 -n/2:-1];
ks=fftshift(k);

u=sech(t);

ut=fft(u);
noise=30;
avg=zeros(1, n);
realizations=30;
for j=1:realizations % 30 sample signals in time and frequency domain, each with periodic domain of T=30
    utn=ut+noise*(randn(1, n)+i*randn(1, n));
    un=ifft(utn);
    avg=avg+utn;
    
    avg2=abs(fftshift(avg))/realizations;
    unf2=ifft(avg);
    subplot(2, 1, 1); plot(t, u, 'k', t, unf2);
    subplot(2, 1, 2); plot(ks, abs(fftshift(ut)), 'r', ks, avg2, 'k');
    
    pause(0.5);
end

unf=ifft(avg);
avg=abs(fftshift(avg))/realizations;

subplot(2, 1, 1); plot(t, u, 'r', t, un, 'c', t, unf, 'k');
subplot(2, 1, 2); plot(ks, abs(fftshift(ut)), 'r', ks, abs(fftshift(utn)), 'c', ks, avg, 'k');

No comments:

Post a Comment