Azzera filtri
Azzera filtri

fractional delay using FFT,IFFT

27 visualizzazioni (ultimi 30 giorni)
zozo
zozo il 22 Nov 2011
Commentato: Shmuel il 11 Mar 2014
Hello, I have been working on delaying any given signal with subsample accuracy (fractional+interger) delay in Frequency domain which results in simple phase change. I know there are functions available in toolboxes(example delayseq), but I would like to do it manually in my program. Here is the code i have written so far:
clc
clear all
close all
Fs=1000;
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T;
delta_T=2.345; %delay time in milliseconds
% Time vector
x = sin(2*pi*50*t);
w=2*pi*50; %angular freq. component
X=fft(x);
Y=X.*exp(-j*w*t*delta_T);
y_1=abs(ifft(Y));
%plot signals
figure;
plot(Fs.*t(1:50),x(1:50))
hold on;
plot(Fs*t(1:50),y_1(1:50),'r');
legend('Original signal','shifted signal');
xlabel('time (milliseconds)')
ylabel('amplitude');
[[Note: I have edited my previous code.]]
My goal is to delay or advance the above signal x (sine) by any amount of time (lets say by 2.345 milliseconds)
Iam getting a weird plot!! :(
Please help me what am I missing?
Thanks!
  10 Commenti
zozo
zozo il 23 Nov 2011
Sir, the flow i wrote above has nothing to do with the code..iam just generalising the steps. My goal is to delay or advance any signal by any amount of time (example: 23.123 ms, 34.742 ms etc.)
shifting the signal by 23ms is easy but the fractional part which is 0.123ms is critical.
Shmuel
Shmuel il 11 Mar 2014
very nice example,
use Y=X.*exp(-j*w*(t+delta_T)); % delay defenition
y_1=imag(ifft(Y)); %% Take image part or real, but not abs.

Accedi per commentare.

Risposta accettata

Walter Roberson
Walter Roberson il 24 Nov 2011
You probably should not be multiplying X(1) by the delay, as X(1) is the DC magnitude that scales the rest of the values.
  3 Commenti
Walter Roberson
Walter Roberson il 25 Nov 2011
I observed some unexpected things that will take me more time to work through. Unfortunately I will probably not be able to work on this until tomorrow now (remote access to my server is still broken!)
Dr. Seis
Dr. Seis il 15 Mag 2012
This answer should be edited or removed. Doing the time-shift correctly (which is not done in the original question) means that the "delay" you say is multiplied to "X(1)" will simply be 1 (the frequency at X(1) is 0, so no "delay" is applied since "X(1)*exp(0)" still equals "X(1)"). Having this as the accepted answer may be misleading.

Accedi per commentare.

Più risposte (2)

Teja Muppirala
Teja Muppirala il 25 Nov 2011
Ok. You've pretty much written exactly what you need to do. This is the right idea:
signal(x)->X=FFT(x)->Y=X*exp(-j*2*pi*f*dt)-->y=ifft(Y)->plot x & y.
But you are multiplying by the wrong factor in your code.
Consider a concrete example. Say t = 0:0.1:0.9 and x = some x(t)
There are 10 values in x, and the period is one.
The Fourier transform gives you the coefficients of the C given below. x(t) = C0 + C1*exp(2i*pi*1*t) + C2*exp(2i*pi*2*t) + ... + C5*exp(2i*pi*5*t) + C(-4)*exp(2i*pi*(-4)*t) + C(-3)*exp(2i*pi*(-3)*t) + ... + C(-1)*exp(2i*pi*(-1)*t)
If instead of x(t) I had x(t+dt), then substituting t --> (t+dt) gives:
x(t+dt) = C0 + C1*exp(2i*pi*1*(t+dt)) + C2*exp(2i*pi*2*(t+dt)) + ... C(-1)*exp(2i*pi*(-1)*t)
= C0 + C1*exp(2i*pi*1*dt)*exp(2*pi*1*t) + C2*exp(2i*pi*2*dt)*exp(2*pi*2*t) + ... C(-1)*exp(2i*pi*(-1)*dt)*exp(2i*pi*(-1)*t)
So you see the coefficients get changed by C = C .* exp(2*pi*[(0:5) (-4:1)]i*dt);
Putting it all together, delay in the time domain is achieved by the appropriate multiplication in the frequency domain:
dt = -0.0521; %<-- Shift by this amount
t = 0:0.01:0.99;
x1 = sin(2*pi*t) + cos(4*pi*t) + sin(6*pi*t);
x2 = sin(2*pi*(t+dt)) + cos(4*pi*(t+dt)) + sin(6*pi*(t+dt));
X1 = fft(x1);
X2 = fft(x2);
X1_DELAY = X1 .* exp(2i*pi*[0:50 -49:-1]*dt);
plot(t,x1,t,x2,t,real(ifft(X1_DELAY)),'r:')
legend({'x(t)' 'x(t+dt)', 'x(t+dt) by FFT'})
norm(X2-X1_DELAY) %<--- Very small
Or to show it using your code:
clear all
close all
Fs=1000;
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T;
delta_T=2.345 / 1000; %delay time in SECONDS
% Time vector
x = sin(2*pi*50*t);
w=2*pi*50; %angular freq. component
X=fft(x);
Y=X.*exp(1j *2*pi*([0:L/2 -L/2+1:-1])*L*T*delta_T);
y_1=real(ifft(Y));
%plot signals
figure;
plot(Fs.*t(1:50),x(1:50))
hold on;
plot(Fs*t(1:50),y_1(1:50),'r');
legend('Original signal','shifted signal');
xlabel('time (milliseconds)')
ylabel('amplitude');
  4 Commenti
michael scheinfeild
michael scheinfeild il 27 Ago 2012
Modificato: michael scheinfeild il 27 Ago 2012
very good example of fft pairs and delay thank you just mention that after fft for plot you can use fft shift that shift the result around the center
i found some issue in the frequency vector if fs is different lets say 1500 and L= 1100 , the function is not correct !! i have also corrected the plot *!!!, should multipy by 1000 for ms clear all close all Fs=6500; T = 1/Fs; % Sample time L = 22000; % Length of signal t = (0:L-1)*T; delta_T= 1 / 1000; %delay time in SECONDS % Time vector fa=250; x = sin(2*pi*fa*t); w=2*pi*fa; %angular freq. component X=fft(x); Y=X.*exp(1j *2*pi([0:L/2 -L/2+1:-1])*L*T*delta_T); y_1=real(ifft(Y)); [sigWithDelay] = delayTheSignal(x,Fs,delta_T); %plot signals figure; plot(1000*t(1:100),x(1:100)) hold on; plot(1000*t(1:100),y_1(1:100),'r'); plot(1000*t(1:100),sigWithDelay(1:100),'ko'); legend('Original signal','shifted signal'); xlabel('time (milliseconds)') ylabel('amplitude'); %================== % signal(x)->X=FFT(x)->Y=X*exp(-j*2*pi*f*dt)-->y=ifft(Y)->
function [sigWithDelay] = delayTheSignal(x,fs,delta_T)
% delta_T = delay seconds
% fs sampling frequency
nfft = 2.^ceil(nextpow2(length(x)));
xfft = fft(x,nfft);
T =1/fs;
resFFT =fs/nfft;
vf = [0:nfft/2 (-nfft/2+1):-1]*resFFT;
Y=xfft.*exp(1j *2*pi*(vf*delta_T));
sigWithDelay =real(ifft(Y));
==========================
Alex
Alex il 1 Nov 2012
Unfortunately, the demo using: x1 = sin(2*pi*t) + cos(4*pi*t) + sin(6*pi*t); is very nice for demo but conceals practical issues. Using demo with signals that wrap nicely gives very good results. When the frequency is changed from 1,2 and 3 to a non integer value (e.g. x1 = sin(2.1*pi*t) + cos(4.3*pi*t) + sin(6.87*pi*t); the outcome will not match any more for a section of time. Now question is why and how to overcome... Here one will need knowledge in signal processing.

Accedi per commentare.


Wayne King
Wayne King il 22 Nov 2011
Not sure if the DSP System Toolbox is an option for you, but if so, please see:
fdesign.fracdelay
  2 Commenti
zozo
zozo il 22 Nov 2011
no sir..i need to implement the algorithm manualy..without using functions..so that i get to learn the actual processing.. thank you.
zozo
zozo il 22 Nov 2011
can anybody help me in this?

Accedi per commentare.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by