Welcome to EDAboard.com

Welcome to our site! EDAboard.com is an international Electronics Discussion Forum focused on EDA software, circuits, schematics, books, theory, papers, asic, pld, 8051, DSP, Network, RF, Analog Design, PCB, Service Manuals... and a whole lot more! To participate you need to register. Registration is free. Click here to register now.

[SOLVED] Implementation of non-traditional FFT

Status
Not open for further replies.

darktronic

Newbie level 6
Joined
Jul 17, 2012
Messages
12
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,368
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: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

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
Joined
Jul 16, 2010
Messages
923
Helped
295
Reputation
590
Reaction score
268
Trophy points
1,343
Activity points
8,543
you haven't defined p(q)
 

darktronic

Newbie level 6
Joined
Jul 17, 2012
Messages
12
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,368
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
Joined
Jul 11, 2012
Messages
281
Helped
48
Reputation
96
Reaction score
47
Trophy points
1,308
Activity points
2,835
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
Joined
Oct 6, 2011
Messages
923
Helped
214
Reputation
428
Reaction score
208
Trophy points
1,333
Location
Bangalore
Activity points
7,438
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
Joined
Jul 17, 2012
Messages
12
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,368
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.
Capture.PNG
 

Mityan

Full Member level 5
Joined
Jul 11, 2012
Messages
281
Helped
48
Reputation
96
Reaction score
47
Trophy points
1,308
Activity points
2,835
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
Joined
Jul 17, 2012
Messages
12
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,368
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(n) with zeroes till it is same length as X(k)?. P(n) in this case is same length as X(k) actually. It is more like x(n).
 

Mityan

Full Member level 5
Joined
Jul 11, 2012
Messages
281
Helped
48
Reputation
96
Reaction score
47
Trophy points
1,308
Activity points
2,835
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
Joined
Jul 17, 2012
Messages
12
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,368
Yes thank you. I understand it now. But p(n) is in this case is not in order. p(n) =randi(N,1,N). So how do I implement it in the FFT?
 

Mityan

Full Member level 5
Joined
Jul 11, 2012
Messages
281
Helped
48
Reputation
96
Reaction score
47
Trophy points
1,308
Activity points
2,835
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(n) or not.
 

darktronic

Newbie level 6
Joined
Jul 17, 2012
Messages
12
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,368
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.

Capture.PNG

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
Joined
Jul 11, 2012
Messages
281
Helped
48
Reputation
96
Reaction score
47
Trophy points
1,308
Activity points
2,835
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
Joined
Jul 17, 2012
Messages
12
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,368
Really? What i'm trying do is based on this formula.

Capture.PNG

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
Joined
Jul 11, 2012
Messages
281
Helped
48
Reputation
96
Reaction score
47
Trophy points
1,308
Activity points
2,835
type (d') instead of d
 

darktronic

Newbie level 6
Joined
Jul 17, 2012
Messages
12
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,368
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
Joined
Jul 11, 2012
Messages
281
Helped
48
Reputation
96
Reaction score
47
Trophy points
1,308
Activity points
2,835
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:2:(n-1)),p(1:2:(n-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
Joined
Jul 17, 2012
Messages
12
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,368
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(n) = exp(-1i*fspacing*p(n)*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
Joined
Jul 11, 2012
Messages
281
Helped
48
Reputation
96
Reaction score
47
Trophy points
1,308
Activity points
2,835
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
Joined
Jul 17, 2012
Messages
12
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,281
Activity points
1,368
Yes, it makes sense actually. If I DFT x(n) 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.

77391d1342608659-capture.png
 

Status
Not open for further replies.

Similar threads

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Top