fractional delay using FFT,IFFT

Question

zozo 2011-11-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/21883-fractional-delay-using-fft-ifft

评论： Shmuel 2014-3-11

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 个评论
显示 8更早的评论隐藏 8更早的评论

zozo 2011-11-23

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 2014-3-11

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.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2011-11-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/21883-fractional-delay-using-fft-ifft#answer_29226

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 个评论
显示 1更早的评论隐藏 1更早的评论

Walter Roberson 2011-11-25

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 2012-5-15

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.

请先登录，再进行评论。

Answer 2

Teja Muppirala 2011-11-25

6
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/21883-fractional-delay-using-fft-ifft#answer_29243

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

michael scheinfeild 2012-8-27

编辑：michael scheinfeild 2012-8-27

在 MATLAB Online 中打开

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 2012-11-1

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.

请先登录，再进行评论。