Prev: position control of x y z direction
Next: Invalid Handle object in GUI
From: Johannes Buechler on 14 Jul 2010 03:36
Rune Allnor <allnor(a)tele.ntnu.no> wrote in message <3dab608e-d1bf-4015-827b-b7fc3f0fbab1(a)r9g2000vbk.googlegroups.com>...

> The factor 2'ish issue (the DC and nyquist components must
> be treated separately, as they only appear once in the spectrum)
> is one of a large number os issues that must be handled correctly
> to get the correct numbers in the end.

correct

> The problem as stated in this tread was how to relate the computed
> numbers to sound pressure. You can't do that without the calibration
> factors.

It would be appreciated a lot if you could finally stop this repeated and unhelpful refering to these ominous "calibration factors"; you have neither at any point clarified what exactly you mean with "calibration data" nor have you tried to help finding the real reasons of the difference between Matlab and other programs.

The problem was solved quite a while ago, the plots are equal now. Since some people might experience similar problems in the future, here the reasons for the (now gone) differences in the Campbell diagrams that this thread was about:

1.) As Godzilla had pointed out somewhere during the thread, the windowing (Hann) needs to be compensated for. I did not do that in my original code (thought I did so though --> thanks Godzilla). Interestingly, the compensation of "2" is obviously applied by both PAK and ArtemiS in all FFT plots, leading to results that show (with broadband noise) levels that are something like 1..2 dB higher than if a rectangular window was used instead. However, the windows correction "2*sqrt(2/3)" is applied in overall level plots. So as to get an overall level plot equal to PAK (ArtemiS uses a different algorithm following those of sound level meters), I have to "re-correct" the sums of each window corrected FFT set by 2*sqrt(2/3)/2, that is, the energy values need to be multiplied by 2/3. I am not sure why PAK uses different windows correction factors depending on the plot to show. For a discussion
of the "correct" correction factor after application of a HANN window, see http://maintenanceforums.com/eve/forums/a/tpc/f/3751089011/m/9711055853?r=3921059853#3921059853

2.) As Mr Käferstein had told me in a PM, the FFT in Matlab needs to be scaled by a factor of 1/(nfft * sqrt(2)) in order to get the expected levels. The advantage of the 1/nfft scaling is that one gets the correct overall level after summing up all the Stützstellen (frequency points) - something that obviously makes perfect sense and is standard in all analysis software that I've see so far. The second part, sqrt(2) is getting the numerical values right if RMS values are to be considered (as usually done since they represent the energy of a signal and are the basis for the SPL calculation).

3.) The setting called "Überlappung" (Overlap) does not mean the same thing in PAK as I had originally assumed. I had expected it to refer to the amount of samples which are considered in both of two successive/adjacent FFT windows, as the window is moved less than its own size. In PAK however, this (FFT block) overlap is only implicitely set via the step size of the reference parameter (time). This step size is, however, (other than the "Überlappung" which is set in the main window of FFT analysis settings) set in a different window which made the whole thing a little confusing to me.

After these three points, numerical values were identical. To get the 3D plot right as well, a last fix had to be done:

4.) The colormap that I used in PAK as approximation of the Matlab "Jet" colormap was incorrect. It had been created visually (=quick&dirty) by a co-worker earlier who had simply compared resulting plots from different programs. Now I have plotted the RGB values in Matlab, resampled them from 64 to 100 (due to a difference in the number of rgb values defining the whole map in both programs) and created a PAK colormap with these values.

The resulting Campbell diagrams now look exactly identical - on my screen up to approx. 2000 Hz. Above 2 kHz the Matlab plots still look different from PAK or ArtemiS even though its numerical values are identical. The latter becomes clear when plotting the averaged FFT, which are identical between all three programs up to fs/2 (or in case of PAK, until where it calculates which is somewhat lower than fs/2). Thus, the difference in the 3D plots could be an issue of different color interpolation for higher frequencies, especially in logarithmic representations of the f-axis. If a solution can be found for that, great - otherwise I can live with that as I know the numerical values are identical (which was my main concern).

Also the overall level can be plotted with identical results between PAK and Matlab.
Here are some plots showing the result. Other than the pictures that were named earlier in this thread, I'll try to keep them available on my domain for a longer time period.

(1) Campbell:
http://www.JohannesBuechler.de/FFT/Campbell_Matlab-vs-PAK.jpg
http://www.JohannesBuechler.de/FFT/Campbell_Matlab-vs-ArtemiS.jpg

(2) Average FFT:
http://www.JohannesBuechler.de/FFT/AverageFFT_Matlab-vs-PAK.jpg

Comparison different HANN window correction factors:
http://www.JohannesBuechler.de/FFT/Matlab_RECT-vs-HANN_CorrectionTypes.jpg

Best regards

Johannes
From: Joerg Braeuer on 22 Jul 2010 10:31
Hello Johannes,

I encountered the same problem when trying to plot Artemis-data in Matlab. If I understood correctly, there are two window-compensation factors. A Division by 2 to compensate for the hann-window (leakage?) and another factor of 2*sqrt(2/3) to get the rms-values, meaning there is an overall scaling factor of sqrt(2/3). Furthermore the scaling factor of 1/(nfft*sqrt(2)) needs to be applied.

If I incorporate these information, my code looks like this:

test_signal = Data{2,1}.channels.values;
Fs = Data{2,1}.Fs;

nfft = 16384;
noverlap = floor(0.5*nfft);
window = hanning(nfft);

[B,F,T] = spectrogram(test_signal,window,noverlap,nfft,Fs);

p0 = 2*10^-5;
X = B*sqrt(2/3)/(nfft*sqrt(2));
L_p = 20*log10(abs(X)/p0)

To plot the campbell, I use:

rpm = Data{3,1}.channels.values;
idx = round((T*Fs+1)); % The indices of the start times and set
rpmT = rpm(idx); % the Abscissa to the speed in RPS -> Campbell
rpmT = (rpmT+[rpmT(2:end),rpmT(end)])/2;

figure('Color',[1 1 1]);
ax1=axes;
surface(rpmT,F,L_p)
shading('interp');
set(gcf,'Renderer','zbuffer');
title('Spectrogram');
xlabel('rpm');
yLabel('frequency [Hz]');
colormap(mycolormap);
view(-90,90); %rotate view
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',[0,350]);
set(ax1,'ydir','reverse');
set(ax1, 'CLim', [0 100]);
colorbar;

However, the matlab-campbell is still way off of what I get using Artemis (about 20 dB). Hopefully, you can help me to understand what I do wrong.

Best regards
Joerg Braeuer
From: Johannes Buechler on 12 Aug 2010 05:45
"Joerg Braeuer" <aballister(a)DELgmx.net> wrote in message <i29kn7$g63$1(a)fred.mathworks.com>...
> Hello Johannes,
>
> I encountered the same problem when trying to plot Artemis-data in Matlab. If I understood correctly, there are two window-compensation factors. A Division by 2 to compensate for the hann-window (leakage?) and another factor of 2*sqrt(2/3) to get the rms-values, meaning there is an overall scaling factor of sqrt(2/3). Furthermore the scaling factor of 1/(nfft*sqrt(2)) needs to be applied.
>
Hi Jörg,

sorry, I only now saw this message, but we're already figured the issue out in PM. For others, here are the necessary code lines to get the right plots:

if strcmpi(Scaling, 'peak')
ScalingFactor = 1;
else %RMS
ScalingFactor = 2^(-1/2);
end

if (strfind(lower(wintype), 'rect'))
window = rectwin(nfft);
WindowCorrectionFactor = 1;
else %HANN
window = hann(nfft);
WindowCorrectionFactor = 2;
end

[B, F, T ] = specgram(signal,nfft,Fs,window,noverlap);

% Scaling
FFTScalingFactor = WindowCorrectionFactor * ScalingFactor / nfft;
B = FFTScalingFactor * abs(B);

B(2:end-1,:) = 2.*B(2:end-1,:); % single side spectrum with even block length

Best regards, Johannes
 | 
Pages: 1
Prev: position control of x y z direction
Next: Invalid Handle object in GUI