lunedì 1 novembre 2010

Time shift of signal using FFT in matlab

Time shift of signal using FFT in matlab:
Let's start with the properties of the fourier trasform:

An integrable function is a function ƒ on the real line that is Lebesgue-measurable and satisfies
\int_{-\infty}^\infty |f(x)| \, dx < \infty.

Basic properties

Given integrable functions f(x), g(x), and h(x) denote their Fourier transforms by \hat{f}(\xi), \hat{g}(\xi), and \hat{h}(\xi) respectively. The Fourier transform has the following basic properties (Pinsky 2002).
Linearity
For any complex numbers a and b, if h(x) = (x) + bg(x), then  \hat{h}(\xi)=a\cdot \hat{f}(\xi) + b\cdot\hat{g}(\xi).
Translation
For any real number x0, if h(x) = ƒ(x − x0), then  \hat{h}(\xi)= e^{-2\pi i x_0\xi }\hat{f}(\xi).
Modulation
For any real number ξ0, if h(x) = e2πixξ0ƒ(x), then  \hat{h}(\xi) = \hat{f}(\xi-\xi_{0}).
Scaling
For a non-zero real number a, if h(x) = ƒ(ax), then  \hat{h}(\xi)=\frac{1}{|a|}\hat{f}\left(\frac{\xi}{a}\right).     The case a = −1 leads to the time-reversal property, which states: if h(x) = ƒ(−x), then  \hat{h}(\xi)=\hat{f}(-\xi).
Conjugation
If h(x)=\overline{f(x)}, then  \hat{h}(\xi) = \overline{\hat{f}(-\xi)}.
In particular, if ƒ is real, then one has the reality condition\hat{f}(-\xi)=\overline{\hat{f}(\xi)}.
And if ƒ is purely imaginary, then  \hat{f}(-\xi)=-\overline{\hat{f}(\xi)}.
Convolution
If h(x)=\left(f*g\right)(x), then   \hat{h}(\xi)=\hat{f}(\xi)\cdot \hat{g}(\xi).


We want to use the traslation property of the fourier trasform to shift a signal:

First of all we want to simply fourier trasform the signal X:

 

TF=fft(X) 

 

Second we have to define the vector of the frequencies used in the exponential that allow the shift:

 

w=[0:floor((numel(X)-1)/2) -ceil((numel(X)-1)/2):-1]/(Xdelta*numel(X));

 

where Xdelta is 1/Fs (Fs is the sample frequency of the signal).

We now define a value of the delay in second (accord to our signal in SI):

 

delay=n;

 

and we apply the sfhit property of the fft: 

 

 TFshifted=TF.*exp(-i*w*2*pi*delay)

 

remember that TF is the fft of the signal X, so a complex function.

We now have the time shifted signal. We want now to build the new signal Xf, that is the correct shifted signal:

 

n1=ceil(abs(delais)/Xdelta);

 

where n1 is the first useful sample in the TFshifted function. This value will be used to build the new signal.

We go back in time domain using the inverse fourier trasform:

 

Xs=real( ifft(TFshifted) );

 

and if delay <0:

 

OUTPUT_SIGNAL(1:numel(X)-n1)=Xs(1:numel(X)-ndeb);

 

else if delay >0:

 

OUTPUT_SIGNAL(ndeb:numel(X))=Xs(ndeb:numel(X));

Here a complete history of the command:

X=signal %signal to shift
Fs=nn %sample frequency
Xdelta=1/Fs
delay=n;
TF=fft(X);
w=[0:floor((numel(X)-1)/2) -ceil((numel(X)-1)/2):-1]/(Xdelta*numel(X));

 TFshifted=TF.*exp(-i*w*2*pi*delay)
n1=ceil(abs(delay)/Xdelta);
Xs=real( ifft(TFshifted) );
if delay <0:
OUTPUT_SIGNAL(1:numel(X)-n1)=Xs(1:numel(X)-n1);
elseif delay >0:
 OUTPUT_SIGNAL(n1:numel(X))=Xs(n1:numel(X));

end

 

HERE A USEFUL FUNCTION THAT YOU CAN DIRECTLY USE:

 

function z=FFTshift(signal,delay,sample_Frequency);
% Created by Piero Poli polipiero85@gmail.com

X=signal; %signal to shift
Fs=sample_Frequency; %sample frequency
Xdelta=1/Fs;

TF=fft(X);
w=[0:floor((numel(X)-1)/2) -ceil((numel(X)-1)/2):-1]/(Xdelta*numel(X));

 TFshifted=TF.*exp(-i*w*2*pi*delay);
n1=ceil(abs(delay)/Xdelta);
Xs=real( ifft(TFshifted) );
if delay <0
OUTPUT_SIGNAL(1:numel(X)-n1)=Xs(1:numel(X)-n1);
elseif delay>0
 OUTPUT_SIGNAL(n1:numel(X))=Xs(n1:numel(X));
end

z=OUTPUT_SIGNAL;

 

the code was tested using a simple dirac function and define the Fs=1. Both positive and negative shift were realized. Some problem are coming out in case of delay similar to the lenght of the signal. The shape of the shifted function is really simila to the input signal, no problem of amplitude are documentated. Some little oscillation can be present in case of really simple function (as dirac pulse), this are due to the fft algorithm in matlab and residual imaginary part from the fourier trasform.

 

If you want to try it use:

 

signal=zeros(1,1000);

signal(500)=100;

z=FFTshift(signal,-100,1); %negative delay

z=FFTshift(signal,100,1); %positive delay

1 commento: