Power spectral density plot converted to sound pressure level
Show older comments
I have analysed a voltage signal and plotted the PSD. However, I need to produce a power-spectral-density plot from the data in terms of non-dimensional frequency. How do I calculate the Strouhal number for the frequency range and generate a power spectral density per unit Strouhal number plot for it?
nfft=2048;
window = hann(nfft);
overlap = 0;
[PSD, freq] = pwelch(det_data, window, overlap, nfft,fs);
figure;
plot(freq, 10*log10(PSD), 'b');
title('Power Spectral Density');
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');
Answers (2)
Hiro Yoshino
on 3 Jun 2023
The simplest way is this; try the following code without receiving the result from the function:
pwelch(det_data, window, overlap, nfft,fs);
It will show the PSD plot with the x-axis set normalized frequency.
Star Strider
on 3 Jun 2023
0 votes
To generate an SPL plot, and assuming compatible units (among other things) —
nfft=2048;
window = hann(nfft);
overlap = 0;
[PSD, freq] = pwelch(det_data, window, overlap, nfft,fs);
figure;
plot(freq, 10*log10(PSD), 'b');
title('Power Spectral Density');
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');
Pref = 20; % microPascals
SPL = pow2db(PSD/Pref^2);
figure
plot(freq, SPL)
grid
Calculating the Strouhal number is not possible.
While we can assume that f is frequency, 
and 
are nowhere defined, miuch less calculated.
and are nowhere defined, miuch less calculated. .
8 Comments
S_saleem
on 3 Jun 2023
Star Strider
on 3 Jun 2023
Yes.
nfft=2048;
window = hann(nfft);
overlap = 0;
[PSD, freq] = pwelch(det_data, window, overlap, nfft,fs);
figure
plot(freq, 10*log10(PSD), 'b')
grid
title('Power Spectral Density')
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density (dB/Hz)')
Pref = 20; % Units: microPascals
SPL = pow2db(PSD/Pref^2);
figure
plot(freq, SPL)
grid
title('Sound Pressure Level')
xlabel('Frequency (Hz)')
ylabel('SPL (Units)')
Uj = 349; % Units: m/s
Dj = 8.08; % Units: mm
St = freq*Dj/Uj; % Strouhal Number
figure
plot(freq, St)
title('Strouhal Number')
xlabel('Frequency (Hz)')
ylabel('Strouhal Number (dimensionless)')
Again, assuming all the units are compatible (currently ‘St’ is in units of ‘cycles/s * mm * s/m = cycles/s * mm * 1E+3m/mm * s/m = cycles*1E+3’ so that may need to be adjusted if it is supposed to be dimensionless) and the pwelch call is correct (specificaly the window size), that should work.
.
S_saleem
on 3 Jun 2023
Please read this "The calculated power spectral density is then converted to sound pressure level (SPL) in decibels (dB) referenced to 20 microPa. The resultant spectra are nondimensionalized to power spectral density per unit Strouhal number"
Now, after obtaining SPL when I try to plot it with St (strouhal number) how do I get dB/St vs St plot?
Secondly, the St values I obtained from the above method are very large it shoulbe of 0.1<St<1.5 range (attaching the figure for reference). What should I do? Please have a look at my code.
I can’t read 'Code..txt' so I’ll have to proceed without it. (I can’t open it onlline or offline.) The error is: foundation::storage::vfs::UnsupportedIRISchemeException
Using my code (above) and adding my own Fourier transform code to it —
Unz = unzip('data..zip')
M1 = readmatrix(Unz{1}, 'FileType','text')
t = M1(:,1) - M1(1,1);
s = M1(:,2);
format long
% tsts = [mean(diff(t)) std(diff(t)) mode(diff(t)) 1/mean(diff(t)) round(1/mean(diff(t)))]
Fsr = round(1/mean(diff(t)))
[sr,tr] = resample(s, t, Fsr); % Resample To A Consistent Sampling Frequency
fs = Fsr;
figure
plot(tr, sr)
grid
xlabel('Time (s)')
ylabel('Amplitude')
Fn = Fsr/2;
L = numel(tr)
NFFT = 2^nextpow2(L);
FTs = fft((sr-mean(sr)).*hann(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, mag2db(abs(FTs(Iv))*2))
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (sB)')
nfft=2048;
window = hann(nfft);
overlap = 0;
[PSD, freq] = pwelch(det_data, window, overlap, nfft,fs);
figure
plot(freq, 10*log10(PSD), 'b')
grid
title('Power Spectral Density')
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density (dB/Hz)')
Pref = 20; % Units: microPascals
SPL = pow2db(PSD/Pref^2);
figure
plot(freq, SPL)
grid
title('Sound Pressure Level')
xlabel('Frequency (Hz)')
ylabel('SPL (Units)')
Uj = 349; % Units: m/s
Dj = 8.08; % Units: mm
St = freq*Dj/Uj; % Strouhal Number
figure
plot(freq, St)
title('Strouhal Number')
xlabel('Frequency (Hz)')
ylabel('Strouhal Number (dimensionless)')
The signal is contaminated with broadband noise, so frequency-selective filtering is not going to help it. One option to deal with it is to use the Savitzky-Golay filter (sgolayfilt). I usually use a third-order polynomial and then change the 'framelen' to get the result I want.
I have no idea what the time vector actually is.
EDIT — (3 Jun 2023 at 23:09)
Just post the file as text rather than attaching it.
.
S_saleem
on 4 Jun 2023
My pleasure!
When I performed the Savitzky-Golay filter on your data (essentially a FIR filter), it was smoothed, however for what you are doing, I doubt that it is impoortant (it is commented-out here). It actually interferes with the PSD calculations
The only important change I made here is to resample the data, since there were significant variations in the sampling intervals in the original data. (I do not understand the reason for the 8 kHz filter, however there is nothing wrong with using it. I changed it to design an IIR filter, since that is more efficient.) The resampling seems to improve the PSD plots over what they were without it, so that may be the only additional processing necessary.
I also do not understand the ‘timestamps’. I subtracted the first value from the rest of them since it makes more sense to me.
Unz = unzip('data..zip');
sample = readmatrix(Unz{1}, 'FileType','text')
timestamps = sample(:, 1);
voltages = sample(:, 2);
Fsr = round(1/mean(diff(timestamps)))
[sr,tr] = resample(voltages, timestamps, Fsr); % Resample To A Consistent Sampling Frequency
fs = Fsr;
voltages = sr;
timestamps = tr - tr(1);
% figure
% subplot(2,1,1)
% plot(timestamps, voltages)
% xlabel('Time (s)')
% ylabel('Amplitude)')
% subplot(2,1,2)
% plot(timestamps, voltages)
% xlim([0.9 1.1])
% xlabel('Time (s)')
% ylabel('Amplitude)')
% sgtitle('Original')
%
% voltages_filtered = sgolayfilt(voltages, 3, 501); % Savitzky-Golay Filtering
%
% figure
% subplot(2,1,1)
% plot(timestamps, voltages_filtered)
% xlabel('Time (s)')
% ylabel('Amplitude)')
% subplot(2,1,2)
% plot(timestamps, voltages_filtered)
% xlim([0.9 1.1])
% xlabel('Time (s)')
% ylabel('Amplitude)')
% sgtitle('Savitzky-Golay Filtered')
% voltages = voltages_filtered;
% fs = 200000; % Hz sample rate
N = length(voltages);
sensitivity=3.6;
pressure= voltages/sensitivity;
max(pressure) % maximum sound pressure level
rms(voltages) % rms of signal
% Low-pass filter
cutoff_frequency = 80000;
filtered_data = lowpass(voltages, cutoff_frequency,fs, 'ImpulseResponse','iir');
% Detrend
mean_removed = detrend(filtered_data , 'constant');
det_data = detrend(mean_removed);
%% lets go with this one
nfft=2048; % number of points
window = hann(nfft); % window vector of length nfft
overlap = 0;
[psd,freq]=pwelch(det_data, window, overlap, nfft,fs);
% Plot the power spectral density
figure
plot(freq,10*log(psd))
title('Power Spectral Density');
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');
% Define the reference pressure
pref = 20e-6;
% Calculate the Sound Pressure Level (SPL) in decibels
spl = 10 * log10(psd/pref^2);
figure;
plot(freq, spl);
title('Sound Pressure Level (SPL)');
xlabel('Frequency (Hz)');
ylabel('SPL (dB)');
% Calculate non-dimensional Frequency and peak frequency
%%% Exit Temperature Solution
rho=1.4; % Specific Heat Ratio
T0=293; % Stagnation temperature
Mj=1.14; % Exit Mach
R=287; % Universal Gas Constant
T_exit = round((T0)*(1+((rho-1)/2)*(Mj^2))^-1); %Exit Temperature
%%% Exit Velocity given Mach number and Exit Temperature Solution
Uj = round(Mj*(rho*R*T_exit)^0.5);
%%% EXIT DIAMETER
D= 8;
Dj= D*1.01;
%% Uj/Dj
ratio= Uj/Dj;
u1=2.40483; % Tam and Tanna (1982)
uc=0.7*Uj ; % convection velocity
ca=346 ; % ambient speed of sound.
b= sqrt(Mj^2-1);
lmda1=(2*u1)/(Dj*b);
L1= (pi*Dj*b)/u1; % wave length
Mc=uc/ca;
theta=120;
cos=cosd(120);
a=(L1*(1-(Mc*cos)));
f= uc/a; % Peak frequency
%%% STROUHAL NUMBER For peak frequency
strouhal_number = (f*Dj)/Uj; %%% 0.65
% %% PSD in terms of NON-DIMENSIONALIZed frequency
St= ((freq*Dj)/(Uj));
plot(St, spl)
title('Sound Pressure Level (SPL)');
xlabel('St');
ylabel('SPL (dB)');
I did not analyse the Strouhal units this time, however I notice that you are incorporating the ‘ca’ value, so the frequency is likely being converted to wavelength, and if the other units are compatible and of the same magnitudes, it will be dimensionless.
.
S_saleem
on 4 Jun 2023
Star Strider
on 4 Jun 2023
My pleasure!
I do not understand what you are doing, particularly with respect to ‘scaling factor’. (I am specifically helping only to get your code to run correctly.)
‘Please let me know how should I scale the ampltiude to a dB/V scaling factor of about 230 for this plot?’
I do not understand what you are doing well enough to respond to that.
.
Categories
Find more on Parametric Spectral Estimation in Help Center and File Exchange
on 4 Jun 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
