Moving signal in frequency domain (MATLAB)

In this post, i am going to show how a moving signal (in time domain) remain unchanged in the frequency domain using MATLAB

The following codes generate a 21 moving signal with the same profile defined by sech(x)

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

slice=0:0.5:10;

[T, S]=meshgrid(t, slice);
[K, S]=meshgrid(ks, slice);

U=sech(T-10*sin(S));
subplot(2, 1, 1);
waterfall(T, S, U), colormap([0 0 0]), view(-15, 70);

[X,Y] = meshgrid(xgv,ygv) replicates the grid vectors xgv and ygv to produce a full grid. This grid is represented by the output coordinate arrays X and Y. The output coordinate arrays X and Y contain copies of the grid vectors xgv and ygv respectively. The sizes of the output arrays are determined by the length of the grid vectors. For grid vectors xgv and ygv of length M and N respectively, X and Y will have N rows and M columns.

U=sech(T-10*sin(S)).*exp(i*0*T) creates the 21 moving signals in time domain in a matrix of size length(slice)*length(t)

waterfall(T, S, U) plots the moving signals as shown in the figure below

The following MATLAB code plot the moving signals in both the time domain and 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);

slice=0:0.5:10;

[T, S]=meshgrid(t, slice);
[K, S]=meshgrid(ks, slice);

U=sech(T-10*sin(S));
subplot(2, 1, 1);
waterfall(T, S, U), colormap([0 0 0]), view(-15, 70);

for j=1:length(slice)
    UT(j, :)=abs(fftshift(fft(U(j, :))));
end


subplot(2, 1, 2);
waterfall(K, S, UT), colormap([0 0 0]), view(-15, 70);

The figure below shows the moving signals in both the time domain and frequency domain, it can be seen that the frequency spectrum of all the moving signals remain unchanged

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

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

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

Denoising via Filtering in MATLAB

In this post, i am going to show how to denoise a signal using a Gaussian filter in MATLAB

The MATLAB code below generates a noisy signal (i.e., un) from the function sech(x) (i.e., u) by adding a white noise via its frequency domain spectrum (ut is the fft of u and utn is the fft of 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=7;
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 figure below shows the resulting plots

The MATLAB code below shows the use of a gaussian filter which is multiplied with the fft of the noisy signal (i.e., utn) to obtained the fft of a filtered noisy signal (i.e., utnf) which is then converted back to the time domain filtered noisy signal (i.e., unf) that is denoised

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

filter=exp(-k.^2); %guassian filter
utnf=filter.*utn;

unf=ifft(utnf);

subplot(2, 1, 1); plot(t, u, 'k', t, un, 'm', t, unf, 'g');
subplot(2, 1, 2); plot(ks, abs(fftshift(ut))/max(abs(fftshift(ut))), 'k', ks, abs(fftshift(utn))/max(abs(fftshift(utn))), 'm', ...
    ks, fftshift(filter), 'b', ...
    ks, abs(fftshift(utnf))/max(abs(fftshift(utnf))), 'g');

axis([-25 25 0 1]);

The figure below shows the resulting plots (where the green line shows the denoised signal in time and frequency domain

Note that in the above MATLAB code it assume that the user knows the frequency they want to look at (in this case it is centered at k=0), therefore the gaussian filter is center at that frequency range (i.e. that guassian is centered at k=0)

The following MATLAB code allows user to vary the width (i.e., sigma) and the location (i.e., the mu) of the Guassian signal

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

sigma=1;
mu=0;

filter=1 / (sigma * (2*pi)^0.5) * exp(-(k-mu).^2/(2*sigma&2)); %guassian filter
utnf=filter.*utn;

unf=ifft(utnf);

subplot(2, 1, 1); plot(t, u, 'k', t, un, 'm', t, unf, 'g');
subplot(2, 1, 2); plot(ks, abs(fftshift(ut))/max(abs(fftshift(ut))), 'k', ks, abs(fftshift(utn))/max(abs(fftshift(utn))), 'm', ...
    ks, fftshift(filter), 'b', ...
    ks, abs(fftshift(utnf))/max(abs(fftshift(utnf))), 'g');

axis([-25 25 0 1]);

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

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

Inverse Fast Fourier Transformation in MATLAB

In this, i am going to show how the inverse fast fourier transformation can be used to find the $n^{th}$ derivative of a hyperbolic secant function sech(x)

Fourier transform has the an amazing property: $\hat{f^{(n)}(k)}=(ik)^n\hat{f(k)}$ where $k$ is the frequency of the Fourier Transform $\hat{f(k)}=\frac{1}{2L} \int_{-L/2}^{L/2} f(x)exp(-i k \pi x/L)dx$ and $f^{(n)}$ is the $n^{th}$ derivative of the function $f(x)$. Therefore $\hat{f^{(n)}(k)}$ is the Fourier transform of the $n^{th}$ derivative of $f(x)$.

Due to this property, we can find the $n^{th}$ derivative of $f(x)$ by first obtaining the $\hat{f^{(n)}(k)}$ then performing an inverse Fourier transform

The sample code below first computes the first and second derivative of sech(x) (i.e. ud and u2d) then compare it against the ones computed using the method described above using inverse Fourier Transform (ifft) (i.e. uds and u2ds) and finally plot ud versus uds

clear all; close all; clc;

L=20;
n=128;

x2=linspace(-L/2, L/2, n+1);

x=x2(1:n);

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

u=sech(x);
ud=-sech(x).*tanh(x);
u2d=sech(x)-2*sech(x).^3;

ut=fft(u);

uds=ifft((i*k).*ut);
u2ds=ifft((i*k).^2.*ut);

ks=fftshift(k);

plot(x, ud, 'r', x, uds, 'mo');

Fast Fourier Transform Sample Code in MATLAB

The following matlab code computes and plots the fourier transform of the gaussian function exp(x^2) (the periodic domain is set to -L / 2 to L / 2) in MATLAB:

clear all; close all; clc;

L=20;
n=128;

x2=linspace(-L/2, L/2, n+1);

x=x2(1:n);

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

u=exp(-x.^2);

ut=fft(u);

plot(fftshift(k), abs(fftshift(ut)));

The MATLAB function fft run in O(n * log(n)) complexity and assumes the periodic domain from 0 to 2*pi, therefore the value vector k is scaled for that the frequency spectrum ut corresponds to (-L/2, L/2). Note that n is power of 2 since the reason fft is able to speed up from O(n^2) to O(n*log(n)) using divide and conquest by dividing the computation into two separate computation each time. linspace() generates (n+1) values evenly from -L/2 to L/2 (inclusive, therefore should be n+1 instead of n, but x takes only the first n values from x2). fftshift rearranges the outputs of fft, fft2, and fftn by moving the zero-frequency component to the center of the array.