FFT is off by 0.4 dB; How can I account for this in my code?

1 view (last 30 days)
I'm trying to generate a plot of a 60s full scale white noise audio clip to compare with plots of the same audio clip generated by two other programs to calibrate the results of my Matlab analysis. At the moment my narrowband plot is showing an average power of 0.4 dB less than the other two plots generated by separate programs, and I cannot for the life of me figure out how to account for that 0.4 dB. I'm getting an average power of -53.37 dB when I should be getting -52.97 dB.
Just for clarity's sake, I do not own the other two programs, another person is running a quick analysis of the same sound file that I'm using for me. I cannot obtain either of these programs and the other person has no knowledge of Matlab coding, so I'm pretty much up a creek here.
The basic overview of my procedure is that I import the noise using wavread, and then perform an FFT using a Hanning window and 50% overlap. Then I use a 10 seconds exponential moving average on each frequency bin to generate a frequency spectrum of the noise file.
I've used three different methods to calculate the FFT: the periodogram function, a simplified fft calculation, and a more complex calculation taking into account the fact that one needn't double the 0 bin or Nyquist Frequency as they only occur once in the FFT. I have also tried the pwelch function, but I am not satisfied with the results that it gives me nor with the restrictions it places on the whole process (for instance, I don't know the method it uses to average all the windowed FFT's together).
Here's my code. If anyone has any insights as to how I may account for the 0.4 dB, it would be greatly appreciated.
clear all
clc
%Import Noise
noise=wavread('Calibration Noise 60s');
%Trim Zeros
idxs=find(noise,1,'first');
idxl=find(noise,1,'last');
noise=noise(idxs:idxl);
%%Hanning Window
nfft=65536;
SR=48000; %Sample Rate samp/sec
w=hann(nfft); %Window
%%Parse
offset=length(w)*.5; % Samples 50% overlap
lfft=floor(length(noise)/offset); %Number of half fft's
limits=ones(1, lfft);
%Start/Stop of Limits
for k=1:length(limits)
limits(k+1)=k*offset;
end
LL=[limits(1), limits(2:end-1)+1];%Start
UL=limits(2:end);%Stop
parse=cell(1, (lfft-2));
for k=1:length(parse)
parse{k}=noise(LL(k):UL(k+1));
end
%%Calculate FFT for each bin
ndft=cell(1,length(parse));
P=cell(1,length(parse));
freq=cell(1,length(parse));
% f=cell(1,length(parse)); %For periodogram analysis
f=0:SR/nfft:SR/2;
fbin=f(2)-f(1);
%FFT Periodogram Function
% for k=1:length(P)
% [P{k},f{k}]=periodogram(parse{k},w,65536,48000);
% end
% f=f{end};
% fbin=f(2)-f(1);
%FFT complex
for k=1:length(P)
N=length(w);
parse{k}=parse{k}.*w;
ndft{k}=fft(parse{k},nfft);
ndft{k}=ndft{k}(1:N/2+1);
P{k}=ndft{k}.*conj(ndft{k})/(nfft.^2);
P{k}(2:end-1)=2*P{k}(2:end-1);
P{k}=P{k}*4/1.1; %account for hanning window
end
%FFT simplified
% for k=1:length(P)
% N=length(w);
% parse{k}=parse{k}.*w;
% ndft{k}=fft(parse{k},nfft);
% ndft{k}=ndft{k}(1:nfft/2+1);
% P{k}=abs(ndft{k})/(N);
% P{k}=P{k}.^2*8/1.1; %account for hanning window
% end
%%Exponential Average
narrowband=P;
T=(nfft/SR);
tau=10;
alpha=1-exp(-.5*T/tau);
narbin=cell(1,length(narrowband{1}));
for k=1:length(narbin)
narbin{k}=ones(1,length(narrowband));
end
for i=1:length(narbin)
for j=1:length(narbin{1})
narbin{i}(j)=narrowband{j}(i);
end
end
narexp=cell(1, length(narbin));
for h=1:length(narbin)
narexp{h}=filter(alpha, [1 alpha-1], narbin{h});
end
psdn=ones(1, length(narexp));
for m=1:length(narexp)
psdn(m)=narexp{m}(end);
end
%%Generate Plot
p=10*log10(psdn);%Narrowband PSD
figure;semilogx(f,p);grid on;axis tight; %Narrowband Plot
xlim([10 30000]);
title(sprintf('PSD of 2013_10_29_09_56_55'),'interpreter','none','FontSize',14);
xlabel('Frequency (Hz)');
ylabel('Power (dB re 1\muPa / Hz)');

Answers (1)

Rick Rosson
Rick Rosson on 11 Oct 2014
Edited: Rick Rosson on 11 Oct 2014

Categories

Find more on Transforms and Spectral Analysis in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!