FFT don't give correct result

Question

Jan-Bert on 1 Jun 2014

0
Link

Direct link to this question

https://nl.mathworks.com/matlabcentral/answers/132021-fft-don-t-give-correct-result

Answered: Jan-Bert on 31 Aug 2014

Hello i'm try to make a few calculations about the the impulse response of a room. One of the calculations is a basic FFT. The measurement i have made with Room EQ Wizard 5. I export their data and load this in matlab. By the impulse response something goes wrong. I don't get the right(and suspected) information.

This is the information i get

This is the information i suspect(printscreen of Audacity)

the link of the Impulse response sample: ImpulseResponse_7.wav

Code loading wavefile with timeplot and frequency plot

%%Wav File
% 2.2675736961451248E-5 // Sample interval (seconds)
dt = 2.2675736961451248E-5;
a = wavread('ImpulseResponse_7.wav');
[m,n] = size(a); t = 0:dt:(m-1)*dt;t=t-1;
figure; set(gcf,'position',[100 225 1000 500]);
plot(t,a);xlim([-.1 1]);
f_s = 1/dt; %nyquist = f_s/2; df = nyquist/(m/2);
%f = -nyquist:df:nyquist;f(f==0) = [];
clear SPL phi f; % correctie 
[SPL,phi,f] = FFT_Amp_Ph(a,f_s,0); 
SPL = SPL(m/2:end); phi = phi(m/2:end); f = f(m/2:end);
figure; set(gcf,'position', [100 400 800 600]);
subplot(2,1,1); semilogx(f,SPL);xlim([10 25000]);
subplot(2,1,2); semilogx(f,phi);xlim([10 25000]);

The code of the modified FFT is:

function [FFT1,FFT2,f] = FFT_Amp_Ph(data,Fs,Spectype)
%Inputs
%   Data: Data array of a time signal 
%   Fs: Sample rate of time signal
%   SpecType: Type of frequency spectrum that return to calc
% Outputs
%   FFT1/ FFT2: 
% Spectype = 0 for Amp/ Phase spectrum; 1 for Real and Imag Spectrum
% Amnplitude/ Phase:
%   FFT1 Returns amplitude of the array (abs from elements of the complex array Z)
%   FFT2 Returns phase of the array (arctan(Im/Re) from elements of the
%   complex array Z)
% Real/Imaginary: 
%   FFT1 Returns real part of the elements of the complex array Z; 
%   FFT2 Returns imaginary part of the elements of the complex array Z; 
FFT = fft(data); f = Fs*linspace(-.5,.5,length(data));
% It is difficult to identify the frequency components by looking at the original signal. 
% Converting to the frequency domain, the discrete Fourier transform of the noisy signal 
% y is found by taking the fast Fourier transform (FFT):
%%Amplitude and Phase
if Spectype == 0
    Amp = abs(FFT);
    Phase = angle(FFT);Phase = Phase/pi*180;
% The angle function can be expressed as angle(z) = imag(log(z)) = atan2(imag(z),real(z)).
    FFT1 = Amp; FFT2 = Phase;
end
%%Real And Imaginary Spectrum
if Spectype ==1
    Re = real(FFT);
    Im = imag(FFT);
% c = complex(a,b) creates a complex output, c, from the two real inputs.
% c = a + bi; c = 2+3i Real part is 2; Imaginary Part is 3i.
    FFT1 = Re; FFT2 = Im;
end

I hope somebody knows what i do wrong or goes wrong

Thanks Jan-Bert

2 Comments
Show NoneHide None

Samuel Hollis on 22 Aug 2014

Hi, have you had any further luck with your findings? I've just started using Room Eq Wizard and wanted to double check its results in matlab.

dpb on 22 Aug 2014

Open in MATLAB Online

The biggest obvious problem is that the frequency vector is screwed up -- sampled at 44.1 kHz so Nyquist frequency is 22+ kHz but the plot has all the energy from 10 kHz up. I didn't try to dig thru the convoluted code but all that's needed for the PSD is

Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
semilogx(f,2*abs(Y(1:NFFT/2+1)))

Sign in to comment.

Sign in to answer this question.

Answer 1

Daniel kiracofe on 23 Aug 2014

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/132021-fft-don-t-give-correct-result#answer_149583

first, "f = Fs*linspace(-.5,.5,length(data));" is not exactly how matlab outputs the frequency vector. To get something like that, you need to use fftshift(fft(data)).

Second, there is an additional scaling that is needed. you can't just do Amp=abs(FFT). the 0 Hz component (i.e. FFT(0)) needs to be divided by n, and the rest of the vector needs to be divided by n/2.

really, matlab's fft() function isn't all that user friendly. For this reason, I never ever call it directly. I use a little wrapper script I wrote, which you can find here (along with some additional details)

http://www.mechanicalvibration.com/Making_matlab_s_fft_functio.html