Wednesday, January 8, 2014

Add noise to the time series signal using Fourier Transform and Frequency Domain Power Spectrum in MATLAB

In this post, i am going to show how to add noise to the time series signal (assume the time series signal is sech(t) in this post) via its frequency domain spectrum in MATLAB

The MATLAB code below plots the time domain (i.e., u) and frequency domain (i.e., ut) of signal sech(t) (using 512 points in time domain data with periodic domain set to T=30)

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);

subplot(2, 1, 1); plot(t, u);
subplot(2, 1, 2); plot(ks, abs(fftshift(ut))/max(abs(fftshift(ut))));

axis([-25 25 0 1]);

The code abs(fftshift(ut))/max(abs(fftshift(ut))) normalize the power spectrum value to be between 0 and 1

The code below shows the resulting plot

Image and video hosting by TinyPic

The code below adds a white noise with normal distribution (noise~N(0, 1)) to the frequency domain (i.e., utn) which is then converted to back to the time domain (i.e. un)

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=1;
utn=ut+noise*(randn(1, n)+i*randn(1, n));
un=ifft(utn);

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

axis([-25 25 0 1]);

The image below shows the resulting plot

Image and video hosting by TinyPic

No comments:

Post a Comment