# [SOLVED]Implementation of non-traditional FFT

Status
Not open for further replies.

#### darktronic

##### Newbie level 6
Hi,

I've been trying to implementing this dft formula using fft:

for w = 1:N
f(w) = 0;
for q = 1:N
f(w) = f(w) + x2(q)*exp(1i*2*pi*p(q)*w/N);
end
end

where f is the frequency response of x2. Signal p is an array with same length as x2.

I found this code where it is a radix-2 fft algorithm to implement a traditional DFT.

function y = fft_rec(x)
n = length(x);
if n == 1
y = x;
else
m = n/2;
y_top = fft_rec(x(1:2n-1)));
y_bottom = fft_rec(x(2:2:n));
d = exp(-2 * pi * i / n) .^ (0:m-1);
z = d .* y_bottom;
y = [ y_top + z , y_top - z ];
end

any ideas how to modified it to make it applicable to my DFT formula above?

Any help will be greatly appreciated. Thanks

#### permute

##### Advanced Member level 3
you haven't defined p(q)

#### darktronic

##### Newbie level 6
Thanks so much for replying

p(q) is a signal with same length as x2(q).

x2(q) = exp(-1i*p(q)*(100)*delay).

So in theory when I DFT x2(q) i would get a peak at frequency delay.

But DFT is so slow so i want to try FFT instead.

- - - Updated - - -

Thanks so much for replying

p(q) is a signal with same length as x2(q).

x2(q) = exp(-1i*p(q)*(100)*delay).

So in theory when I DFT x2(q) i would get a peak at frequency delay.

But DFT is so slow so i want to try FFT instead.

#### Mityan

##### Full Member level 5
For direct furier transform there is a 'dftmtx' operator in MATLAB. Try to read help on this, it may be useful

#### milind.a.kulkarni

##### Advanced Member level 3
One of the key thing that you need to clear that you need to send the signal sample length in Power of 2 for radix-2 FFT algorithm like 2 ,4,8,16,32.....where as in normal DFT you can have any length for signal...

Good Luck

#### darktronic

##### Newbie level 6
Thank you so much for everyone that replied. I dont think dftmx is going be helpful as it is same as FFT function right?. And I know radix-2 only applicable to sample size N = 2^a where a is integer more than 0. Btw, i attached an image to be more clear with my question. I'm trying to implement FFT on equation 2 but i found that radix-2 algorithm usually applicable only to equation 1. Any ideas on how to implement fft on equation 2? Any thought will be greatly appreciated.

#### Mityan

##### Full Member level 5
dftmtx is not the same as fft function. This command creates DFT matrix of any desired size not only power of 2. For implementing a transfor, you should multiply your vector with that matrix and that means you are not using fast algorithm but the DFT as it is.

By the way in your second equation there is no "-" sign, so it seems to be an inverse DFT.

And at last if your pn values are not from the 0,1,...N-1 set, so I guess you should scale your N and interpolate your pn vector by zeroes to get the following:
K*pn are the elements of 0,1,...K*N-1 set.

#### darktronic

##### Newbie level 6
Yes thank you. I agree dftmtx is not same as fft function, but how do I use it to implement equation 2. Btw, can you explain what 'K*pn are the elements of 0,1,...K*N-1 set' means? Did you mean that padded p with zeroes till it is same length as X(k)?. P in this case is same length as X(k) actually. It is more like x.

#### Mityan

##### Full Member level 5
Your equations 1 and 2 are distinguished only in 1 element - pn instead of n in the exponent.
In the input data vector (time domain) all elements are the integer multiples of sampling period T.
After FFT you get your signal in frequency domain, where frequency samples are integer multiples of (1/T).

So (for example) for N=4 n = 0,1,2,3. And if you have pn = 0,0.5,2,3 instead, you should make N=8 and pn will be 0,1,0,0,4,0,6,0.
This means the interpolation of your signal and the new reduced time step will be the GCD of your pn samples. Got it?

#### darktronic

##### Newbie level 6
Yes thank you. I understand it now. But p is in this case is not in order. p =randi(N,1,N). So how do I implement it in the FFT?

#### Mityan

##### Full Member level 5
I'm not sure, but for every Xk sample you have the sum for n from 0 to N-1, so its no matter - in order p or not.

#### darktronic

##### Newbie level 6
Yeah, but I still can't figure out how to do it. . I got this equation based on radix-2 method, but don't have any idea how to implement it recursively.

This is what I've done but it doesn't work. Can someone point me what is wrong is it? Thanks

Code:
function y = fft_rec(x,p)
n = length(x);
if n == 1
y = x;
else
m = n/2;
y_top = fft_rec(x(1:2:(n-1)),p(1:2:(n-1)));
y_bottom = fft_rec(x(2:2:n),p(2:2:n));
d = exp(2 * pi * 1i / n) .^ p(1:m);
z = d .* y_bottom;
y = [ y_top + z , y_top - z ];
end

#### Mityan

##### Full Member level 5
In formula that is posted in a picture I dont see any recursion, because sample Xk doesnt depend on Xk-1.

#### darktronic

##### Newbie level 6
Really? What i'm trying do is based on this formula.

that is successfully implemented using this code

Code:
function y = fft_rec(x)
n = length(x);
if n == 1
y = x;
else
m = n/2;
y_top = fft_rec(x(1:2:(n-1)));
y_bottom = fft_rec(x(2:2:n));
d = exp(-2 * pi * i / n) .^ (0:m-1);
z = d .* y_bottom;
y = [ y_top + z , y_top - z ];
end

I merely want to modify it to implement my formula the post before

Last edited:

#### Mityan

##### Full Member level 5
type (d') instead of d

#### darktronic

##### Newbie level 6
I'm sorry, but why transposing d would help? I tried it and the code can't run after that.

#### Mityan

##### Full Member level 5
Im sorry but in my MATLAB R2006a it works.

p = randint(128,1,5); or p = rand(128,1);
x = rand(128,1);

With the function

function y = fft_rec2(x,p)
n = length(x);
if n == 1
y = x;
else
m = n/2;
y_top = fft_rec2(x(1:2n-1)),p(1:2n-1)));
y_bottom = fft_rec2(x(2:2:n),p(2:2:n));
d = exp(2 * pi * 1i / n) .^ p(1:m);
z = (d') .* y_bottom;
y = [ y_top + z , y_top - z ];
end

it successfully works with no errors.

#### darktronic

##### Newbie level 6
Oh I got what you mean, you made x and p as row vectors but actually my x and p is column vectors. My code works without error too but it's not what I expected. x actually is x = exp(-1i*fspacing*p*delay). So whenever I DFT it, I would get a peak at f = delay. But my code above doesnt produce the peak so I know it doesn't work actually

#### Mityan

##### Full Member level 5
You sample your signal at a variable random rate. Why you decided that you should have this exact peak. It seems to me sensless

#### darktronic

##### Newbie level 6
Yes, it makes sense actually. If I DFT x using equation 2 below I would get a peak at f=delay. Right now I'm trying to figure out the FFT for that equation. Do you think it's possible?. I've read that radix-2 utilised twiddle factors. But I can't get any twiddle factor out of the equation 2 to implement FFT. So I'm thinking this is probably impossible.

Status
Not open for further replies.