fractional delay using FFT,IFFT
27 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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
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.
Risposta accettata
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
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
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.
Più risposte (2)
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
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
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.
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
Vedere anche
Categorie
Scopri di più su Signal Operations in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!