FFT don't give correct result

6 views (last 30 days)
Jan-Bert
Jan-Bert on 1 Jun 2014
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
dpb
dpb on 22 Aug 2014
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.

Accepted Answer

Daniel kiracofe
Daniel kiracofe on 23 Aug 2014
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)

More Answers (1)

Jan-Bert
Jan-Bert on 31 Aug 2014
Thanks that was the problem.!! It works now fine.
Thanks for the script i made a pdf print of it ;)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by