Incorrect Sound Pressure Level from FFT in Matlab
in [Matlab]
|
From: Johannes Buechler on 8 Apr 2010 09:47 Hello folks, Obviously I am doing something wrong but I can't figure out what it is. My procedure: Having time data measured with a commercial data aquisition tool (in this case, PAK from Muller-BBM). Now I would like to draw a correct Campbell diagram from that. It works and looks more or less exactly like in the professional software, however I have wrong numerical dB values. That is, I have to shift the color axis in order to get the same output. The values in Matlab are incorrect, whereas the professional software shows the correct scaling. The same happens if I take just one block and calculate its FFT and compare, or if I calculate the Overall Level from each spetrum. Does anyone know what could be the reason? The Output for comparison is to be seen here, beginning one hour after this post (because I'll have to upload the files at home, at work it's blocked) http://www.JohannesBuechler.de/Campbell_PAKvsMatlab_v2.bmp http://www.JohannesBuechler.de/Campbell_PAKvsMatlab_v3.bmp v2: SPL calculated with the correct reference value for Sound Pressure Levels of po = 2*10^-5 Pa. It has approximately 50 dB shift relative to the correct graph from PAK. v3: SPL calculated with a reference value of 1 Pa instead of po = 2*10^-5 Pa. This one has not as much shift as v3, but still approximately 20 dB shift. What the heck am I doing wrong? Here's my code: OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO function F = JAB_Campbell2(Y) nfft = 2^14; Fs = 2^14; OvlpPerc = 50; rpm = Y(:,end); signal = Y(:,1); window = hanning(nfft); noverlap = floor(OvlpPerc/100*nfft); % noverlap = 0; effective_BL = nfft - noverlap; NumberOfBlocks = floor(length(signal)/effective_BL); LengthDifference = NumberOfBlocks * effective_BL - length(signal); if (~exist('Weighting','var')) Weighting = 'lin'; end dBText = ['dB(' Weighting ') SPL']; % DEBUG % disp([' Länge der Daten: ' num2str(length(signal)) ' samples']); % % disp(['=> Anzahl Blöcke ohne Overlap wäre: ' num2str(floor(length(signal)/nfft))]); % disp([' Anzahl Blöcke mit ' num2str(OvlpPerc) '% Overlap: ' num2str(NumberOfBlocks) ]); % % disp([' Das entspricht ' num2str(NumberOfBlocks * effective_BL) ' samples']); % disp(['=> Längendifferenz: ' num2str(-LengthDifference) ]); % DEBUG Ende [B,F,T] = specgram(signal,nfft,Fs,window,noverlap); % B ist das Einseitenspektrum inkl. f=0 % => F(1) = 0 % % T enthält die Zeitpunkte der Anfänge jedes Zeitblocks % => T(1) = 0 if (max(size(B))-1 ~= nfft/2) error('Length of calculated FFT vectors incorrect!'); end idx = round(T*Fs+1); % start indices of the (possibly overlapping) FFT blocks for i=1:length(idx) rpmT(i) = mean(rpm(idx(i):idx(i)+nfft-1)); end p0 = 2*10^-5; % L_p = 20*log10(abs(B)/(sqrt(2)*p0)); % v1 L_p = 20*log10(abs(B)/p0); % v2 % L_p = 20*log10(abs(B)); % v3 % disp(' -------------- DEBUG '); % min(min(abs(B))) % max(max(abs(B))) % min(min(L_p)) % max(max(L_p)) % disp(' -------------- DEBUG '); fig1 = figure('Color',[1 1 1]); ax1=axes('XGrid','on','YGrid','on'); surface(rpmT,F,L_p); shading('interp'); set(gcf,'Renderer','zbuffer'); title('Campbell Diagram'); xlabel('Engine Revolution RPM'); ylabel('Frequency [Hz]'); zlabel('TEST'); axis tight; set(ax1,'Layer','top'); set(ax1,'YDir','normal'); set(ax1,'Box','on'); set(ax1,'XMinorTick','on'); set(ax1,'YMinorTick','on'); % set(ax1,'YLim',[20,Fs/2]); set(ax1, 'YLim',[20,2000]); set(ax1, 'YScale', 'log'); % set(ax1, 'XLim', [rpmT(1),rpmT(end)]); set(ax1, 'XLim', [1000,4000]); % set(ax1, 'CLim', [round(min(min(L_p))/5)*5 round(max(max(L_p))/5)*5]); set(ax1, 'CLim', [25 round(max(max(L_p))/5)*5]); % set(ax1, 'CLim', [25 75]); % set(ax1, 'ZLim', [50 150]); % colormap('jet'); c_ax1 = colorbar( ... 'Box','on',... 'YLim',[0 150],... 'Location','SouthOutside' ... ); %% Create textbox annotation1 = annotation(... fig1,'textbox',... 'Position',[0.4393 0.045 0.1607 0.0746],... 'LineStyle','none',... 'String',dBText,... 'FitHeightToText','on'); for i=1:min(size(B)) % für jedes Einzelspektrum SumE(i) = sum(B(2:end,i).^2)/(nfft/2); end MeanE = SumE ./ T(2); % T(2) is the time length of each block % fig2 = figure; % plot(rpmT,20*log10(abs(MeanE))); OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO Thanks for your help :-) Best regards, Joe
From: Johannes Buechler on 8 Apr 2010 11:17 pictures are uploaded now... "Johannes Buechler" <jb_spam(a)gmx.DELnet> wrote in message <hpkmpe$bqf$1(a)fred.mathworks.com>... > Hello folks, > > Obviously I am doing something wrong but I can't figure out what it is. > > My procedure: Having time data measured with a commercial data aquisition tool (in this case, PAK from Muller-BBM). > > Now I would like to draw a correct Campbell diagram from that. It works and looks more or less exactly like in the professional software, however I have wrong numerical dB values. That is, I have to shift the color axis in order to get the same output. The values in Matlab are incorrect, whereas the professional software shows the correct scaling. > > The same happens if I take just one block and calculate its FFT and compare, or if I calculate the Overall Level from each spetrum. > > Does anyone know what could be the reason? > > The Output for comparison is to be seen here, beginning one hour after this post (because I'll have to upload the files at home, at work it's blocked) > http://www.JohannesBuechler.de/Campbell_PAKvsMatlab_v2.bmp > http://www.JohannesBuechler.de/Campbell_PAKvsMatlab_v3.bmp > > v2: SPL calculated with the correct reference value for Sound Pressure Levels of po = 2*10^-5 Pa. It has approximately 50 dB shift relative to the correct graph from PAK. > > v3: SPL calculated with a reference value of 1 Pa instead of po = 2*10^-5 Pa. This one has not as much shift as v3, but still approximately 20 dB shift. > > What the heck am I doing wrong? > > Here's my code: > OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO > function F = JAB_Campbell2(Y) > > nfft = 2^14; > Fs = 2^14; > OvlpPerc = 50; > > rpm = Y(:,end); > signal = Y(:,1); > > window = hanning(nfft); > noverlap = floor(OvlpPerc/100*nfft); > % noverlap = 0; > > effective_BL = nfft - noverlap; > NumberOfBlocks = floor(length(signal)/effective_BL); > LengthDifference = NumberOfBlocks * effective_BL - length(signal); > > if (~exist('Weighting','var')) > Weighting = 'lin'; > end > > dBText = ['dB(' Weighting ') SPL']; > > % DEBUG > % disp([' Länge der Daten: ' num2str(length(signal)) ' samples']); > % > % disp(['=> Anzahl Blöcke ohne Overlap wäre: ' num2str(floor(length(signal)/nfft))]); > % disp([' Anzahl Blöcke mit ' num2str(OvlpPerc) '% Overlap: ' num2str(NumberOfBlocks) ]); > % > % disp([' Das entspricht ' num2str(NumberOfBlocks * effective_BL) ' samples']); > % disp(['=> Längendifferenz: ' num2str(-LengthDifference) ]); > % DEBUG Ende > > [B,F,T] = specgram(signal,nfft,Fs,window,noverlap); > > % B ist das Einseitenspektrum inkl. f=0 > % => F(1) = 0 > % > % T enthält die Zeitpunkte der Anfänge jedes Zeitblocks > % => T(1) = 0 > > if (max(size(B))-1 ~= nfft/2) > error('Length of calculated FFT vectors incorrect!'); > end > > idx = round(T*Fs+1); % start indices of the (possibly overlapping) FFT blocks > > for i=1:length(idx) > rpmT(i) = mean(rpm(idx(i):idx(i)+nfft-1)); > end > > p0 = 2*10^-5; > % L_p = 20*log10(abs(B)/(sqrt(2)*p0)); % v1 > L_p = 20*log10(abs(B)/p0); % v2 > % L_p = 20*log10(abs(B)); % v3 > > % disp(' -------------- DEBUG '); > % min(min(abs(B))) > % max(max(abs(B))) > % min(min(L_p)) > % max(max(L_p)) > % disp(' -------------- DEBUG '); > > fig1 = figure('Color',[1 1 1]); > ax1=axes('XGrid','on','YGrid','on'); > surface(rpmT,F,L_p); > shading('interp'); > set(gcf,'Renderer','zbuffer'); > title('Campbell Diagram'); > xlabel('Engine Revolution RPM'); > ylabel('Frequency [Hz]'); > zlabel('TEST'); > axis tight; > set(ax1,'Layer','top'); > set(ax1,'YDir','normal'); > set(ax1,'Box','on'); > set(ax1,'XMinorTick','on'); > set(ax1,'YMinorTick','on'); > % set(ax1,'YLim',[20,Fs/2]); > set(ax1, 'YLim',[20,2000]); > set(ax1, 'YScale', 'log'); > % set(ax1, 'XLim', [rpmT(1),rpmT(end)]); > set(ax1, 'XLim', [1000,4000]); > % set(ax1, 'CLim', [round(min(min(L_p))/5)*5 round(max(max(L_p))/5)*5]); > set(ax1, 'CLim', [25 round(max(max(L_p))/5)*5]); > % set(ax1, 'CLim', [25 75]); > % set(ax1, 'ZLim', [50 150]); > > % colormap('jet'); > > c_ax1 = colorbar( ... > 'Box','on',... > 'YLim',[0 150],... > 'Location','SouthOutside' ... > ); > > %% Create textbox > annotation1 = annotation(... > fig1,'textbox',... > 'Position',[0.4393 0.045 0.1607 0.0746],... > 'LineStyle','none',... > 'String',dBText,... > 'FitHeightToText','on'); > > > for i=1:min(size(B)) % für jedes Einzelspektrum > SumE(i) = sum(B(2:end,i).^2)/(nfft/2); > end > > MeanE = SumE ./ T(2); % T(2) is the time length of each block > > % fig2 = figure; > % plot(rpmT,20*log10(abs(MeanE))); > > OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO > > Thanks for your help :-) > > Best regards, Joe
From: Rune Allnor on 8 Apr 2010 11:46 On 8 apr, 15:47, "Johannes Buechler" <jb_s...(a)gmx.DELnet> wrote: > Hello folks, > > Obviously I am doing something wrong but I can't figure out what it is. > > My procedure: Having time data measured with a commercial data aquisition tool (in this case, PAK from Muller-BBM). > > Now I would like to draw a correct Campbell diagram from that. It works and looks more or less exactly like in the professional software, however I have wrong numerical dB values. If you want the correct numbers you need to calibrate your measurement system. That's a not at all trivial task, which requires you to determine the (possibly frequency dependent and nonlinear) scaling coefficients throughout your measurement system. Rune
From: Johannes Buechler on 8 Apr 2010 12:03 Rune Allnor <allnor(a)tele.ntnu.no> wrote in message <06b4754e-1c5d-4985- > If you want the correct numbers you need to calibrate > your measurement system. That's a not at all trivial > task, which requires you to determine the (possibly frequency > dependent and nonlinear) scaling coefficients throughout > your measurement system. > > Rune Hi Rune, the system is calibrated of course. This is why the graph from the measurement system shows the correct pressure levels. The same values are exported, so they should also be correct in the Matlab plot. The question is, what is wrong in the Matlab calculation? Anyway, thanks for the hint. Joe
From: Mark Shore on 8 Apr 2010 12:06 Rune Allnor <allnor(a)tele.ntnu.no> wrote in message <06b4754e-1c5d-4985-a4b1-e3416edf7e19(a)r1g2000yqb.googlegroups.com>... > On 8 apr, 15:47, "Johannes Buechler" <jb_s...(a)gmx.DELnet> wrote: > > Hello folks, > > > > Obviously I am doing something wrong but I can't figure out what it is. > > > > My procedure: Having time data measured with a commercial data aquisition tool (in this case, PAK from Muller-BBM). > > > > Now I would like to draw a correct Campbell diagram from that. It works and looks more or less exactly like in the professional software, however I have wrong numerical dB values. > > If you want the correct numbers you need to calibrate > your measurement system. That's a not at all trivial > task, which requires you to determine the (possibly frequency > dependent and nonlinear) scaling coefficients throughout > your measurement system. > > Rune As Rune notes, it appears your code derives dB values through intermediate steps from Y -> signal -> B -> L_p without any explicit calibration constant used. Luckily - based on the very similar appearance of the paired plot - the scaling factor appears to be constant over the range of frequencies you're measuring. Presumably you can derive this by either looking through instrument documentation, or normalizing your MATLAB results to the PAK ones.
|
Next
|
Last
Pages: 1 2 3 4 5 6 Prev: fig cannot open, Error using ==> open at 164 Next: Static text with horizontal slider |