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] Problem with Efcient computation of the DFT of a 2N-point real sequence

Status
Not open for further replies.

karanbanthia

Junior Member level 1
Joined
Jan 7, 2010
Messages
16
Helped
1
Reputation
2
Reaction score
1
Trophy points
1,283
Location
pune
Activity points
1,412
Hi all,
I am trying to compute 2N-point DFT of 2Z-point real valued input signal from N-point FFT, but I am getting mirror image in the first half of the spectrum as well.
Below is my MATLAB code
Input signal consists of frequencies: 40, 400 and 800Hz. Sampling freq. is 3072Hz.
Method used is mentioned in Texas instruments application report 'SPRA291', page no. 9

Code:
% This script computers 2N point FFT from 2N point real input signal by
% generating a complex signal of length N and computing N point FFT of this
% complex signal

Fs = 3072;   % Sampling frequency in Hz
T = 1/Fs;    % Sampling time interval
N = 1024;     % No of samples to be taken
t = (0:T:(N*T)-T);                           % Time vector

% Generate the 2N point real input vibration signal
vibrationSignal = 0.4*sin(2*pi*40*t) + 0.2*sin(2*pi*400*t) + 0.25*sin(2*pi*800*t);
signalX = zeros(1,(N/2));
signalY = zeros(1,(N/2));
complexSignal = zeros(1,(N/2));
finalSpectrum = zeros(1,(N));
A = zeros(1,(N/2));
B = zeros(1,(N/2));
vibration_spectrum_real = zeros(1,(N/2));
vibration_spectrum_imag = zeros(1,(N/2));
k = 1;
for k=1:1:(N/2)
    A(k) = (1 - (1i*exp((-1i*2*pi*k)/N)))/2;
    B(k) = (1 + (1i*exp((-1i*2*pi*k)/N)))/2;
end
% Generate N point signal from 2N point real input signal
n = 1;
for i=1:2:(N-1)
    signalX(n) = vibrationSignal(i);
    signalY(n) = vibrationSignal(i + 1);
    n = n + 1;
end
% Generate N point complex signal
for k=1:1:(N/2)
    complexSignal(k) = signalX(k) + 1i*signalY(k);
end
% Compute FFT of N point complex signal
spectrum1_output = fft(complexSignal,(N/2))/(N/2);
% Perform the split operation
k = 1;
for k=1:1:((N/2))
    vibration_spectrum_real(k) = real(spectrum1_output(k))*real(A(k)) - ...
        imag(spectrum1_output(k))*imag(A(k)) + real(spectrum1_output((N/2)-k+1))*real(B(k)) +...
        imag(spectrum1_output((N/2)-k+1))*imag(B(k));
    vibration_spectrum_imag(k) = imag(spectrum1_output(k))*real(A(k)) + ...
        real(spectrum1_output(k))*imag(A(k)) + real(spectrum1_output((N/2)-k+1))*imag(B(k)) -...
        imag(spectrum1_output((N/2)-k+1))*real(B(k));
    
    vibration_spectrum_real(N-k+1) = vibration_spectrum_real(k);
    vibration_spectrum_imag(N-k+1) = -vibration_spectrum_imag(k);

end

for k=1:1:((N))
    finalSpectrum(k) = vibration_spectrum_real(k) + 1i*vibration_spectrum_imag(k);
end

finalSpectrum = abs(finalSpectrum);
% Convert the x-axis range into Hz
x = ((Fs/N)*0) :(Fs/N) : (Fs/2)-1;
% Plot only the first half of spectrum
plot(x',finalSpectrum(1:N/2));

2N-point real FFT.JPG

As can be seen from the output image, I am getting mirror image of input frequencies around 1536Hz which is unwanted.
Has anyone used this method to calculate 2N point FFT?
 

KlausST

Super Moderator
Staff member
Joined
Apr 17, 2014
Messages
19,455
Helped
4,300
Reputation
8,605
Reaction score
4,260
Trophy points
1,393
Activity points
129,016
Hi,

as far as i can see you are not generating a sinewave that fits (with integer numbers) into the 1024 points.

a 40Hz signal needs 25ms to repeat.

1024/3072 = 333.33ms this is not an integer multiple of 25ms.

You need to adjust sampling frequency, samples count and signal frequencies.

Klaus

Klaus
 

karanbanthia

Junior Member level 1
Joined
Jan 7, 2010
Messages
16
Helped
1
Reputation
2
Reaction score
1
Trophy points
1,283
Location
pune
Activity points
1,412
Hi Klaus,
My input frequency will keep on varying, while sampling freq and number of samples will be const. As per theroy, FFT needs minimum number of samples of each freq. which is satisfied since all the cycle are repeated multiple time in the acquired time frame, so that should not be an issue.
In addition, 512-point complex FFT of the same signal gives perfect results, while 1024-point FFT of the same signal calculated as shown, gives mirror freq. around Fs/4, which is odd. Mirror freq. around Fs/2 is obvious.

Regards,
Karan
 

KlausST

Super Moderator
Staff member
Joined
Apr 17, 2014
Messages
19,455
Helped
4,300
Reputation
8,605
Reaction score
4,260
Trophy points
1,393
Activity points
129,016
Hi,

which is satisfied since all the cycle are repeated multiple time in the acquired time frame
i can´t relate to this

if you have non integer multiples you usually need a window function.

***
for an FFT you need 2^n samples. while for a DFT you you may use any number of samples.

Klaus
 

karanbanthia

Junior Member level 1
Joined
Jan 7, 2010
Messages
16
Helped
1
Reputation
2
Reaction score
1
Trophy points
1,283
Location
pune
Activity points
1,412
Hi Klaus,
Tried windowing as well. It only reduces the spectrum leakage (what it is meant to do). In case no windowing is performed, even though the input signal is non-integer multiples, leakage will be more but it wont cause aliasing around Fs/4.
In my results, aliasing is present around Fs/4 in addition to Fs/2.

Regards,
Karan
 

FvM

Super Moderator
Staff member
Joined
Jan 22, 2008
Messages
48,294
Helped
14,229
Reputation
28,719
Reaction score
12,923
Trophy points
1,393
Location
Bochum, Germany
Activity points
279,648
In case no windowing is performed, even though the input signal is non-integer multiples, leakage will be more but it wont cause aliasing around Fs/4.
Yes. There most be something wrong with the implementation so that the conjugate symmetry doesn't work.
 

FvM

Super Moderator
Staff member
Joined
Jan 22, 2008
Messages
48,294
Helped
14,229
Reputation
28,719
Reaction score
12,923
Trophy points
1,393
Location
Bochum, Germany
Activity points
279,648
To help others who may stumble upon the same problem, can you briefly describe where the error in the original code is?
 

Status
Not open for further replies.

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Top