FFT don't give correct result
11 次查看(过去 30 天)
显示 更早的评论
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/163643/image.png)
This is the information i suspect(printscreen of Audacity)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/163644/image.png)
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 个评论
Samuel Hollis
2014-8-22
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
2014-8-22
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)))
采纳的回答
Daniel kiracofe
2014-8-23
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)
0 个评论
更多回答(1 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!