[SOLVED] Problem with Efcient computation of the DFT of a 2N-point real sequence

[karanbanthia]

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

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

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]

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]

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]

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]

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]

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.

 

[karanbanthia]

Here is the link which helped me solve the problem:
https://simon.rozman.si/computers/dsp/dftreal

Regards,
Karan M. Banthia

 

[FvM]

To help others who may stumble upon the same problem, can you briefly describe where the error in the original code is?