fractional delay using FFT,IFFT

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!

[Walter Roberson]

What is this expression intended to mean:

[0:((numel(X)-1)/2) -((numel(X)-1)/2):-1]

Looks to me to be equivalent to 0:0:-1 which in term would be the empty list because the first element is greater than the last and the increment is not negative.

[zozo]

it is:

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

[zozo]

please help experts!

[Walter Roberson]

Dang spacing in that w expression. It would be clearer if you put a comma between the two parts!

Your code does not work as posted even with the delay of -2. n1 comes out as 20000 but numel(X) is 101, and 1:(101-20000) is the empty array so your OUTPUT_SIGNAL comes out empty if it did not already exist, and unchanged if it already existed.

The delay > 0 case has similar problems. And you have no case for delay of exactly 0, which would leave OUTPUT_SIGNAL undefined.

[zozo]

@Walter: ok..i got your point..after setting the break point in my code, I noticed that my code does not work.!! :(

[zozo]

Can you please provide the theoretical steps required to produce delayed versions of any signal (both interger and fractional multiples of samples).

weblinks/literature is also welcome.

thank you.

[zozo]

signal(x)->X=FFT(x)->Y=X*exp(-j*2*pi*f*dt)-->y=ifft(Y)->plot x & y.

where,

dt=time delay(which can take both integer and fractional values)

f=freq. of signal

Am I missing something?

[Walter Roberson]

Your flow there doesn't appear to have reference to OUTPUT_SIGNAL so I do not know what it is intended to represent.

When you are calculating the shifted fft, you are using delay as number of samples (I think), but when you calculate n1, you are using delay as number of seconds.

[zozo]

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]

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]

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.

[zozo]

so can u plz suggest the change in code?

[Walter Roberson]

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]

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]

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

[zozo]

@Teja: thank you so much for ur clear explaination!! :)

[zozo]

@Teja: If I do the above process efficiently, by performing 'N' point FFT where 'N' is a power of 2, then what changes should be made? (There is an error of order 10^-3 between x2 and X1_DELAY when I replace 'L' with 'N' in the above process.

[michael scheinfeild]

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]

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.

Answer 3

[Wayne King]

Not sure if the DSP System Toolbox is an option for you, but if so, please see:

fdesign.fracdelay

[zozo]

no sir..i need to implement the algorithm manualy..without using functions..so that i get to learn the actual processing.. thank you.

[zozo]

can anybody help me in this?