Power spectral density plot converted to sound pressure level

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)');

Answer 1

[Hiro Yoshino]

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.

[S_saleem]

I need to convert this PSD into SPL and then plot dB/St vs St plot. St=strouhal number. Can you help me with it?

Answer 2

[Star Strider]

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.

.

[S_saleem]

Uj= 349 ms-1 and Dj= 8.08 mm. Is it possible to calculate strohual number now with this?

St= ((freq*Dj)/Uj);

[Star Strider]

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]

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.

[Star Strider]

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')

Unz = 1×1 cell array
    {'sample.lvm'}

M1 = readmatrix(Unz{1}, 'FileType','text')

M1 = 500000×2

1.0e+09 * 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000

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)))

Fsr = 
      200000

[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)

L = 
      499999

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);

Unrecognized function or variable 'det_data'.

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]

I really appreciate your help. Here is my code. The data (sample 500,000 points) is provided in terms of Timestamps vs Voltages. How to perform Savitzky-Golay filter on this data? Is it necessary to apply this filter before analysing power spectral density for it?

sample= load('sample');
timestamps = sample(:, 1);
voltages = sample(:, 2);
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);

% 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
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)');

[Star Strider]

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')

sample = 500000×2

1.0e+09 * 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000 3.6734 0.0000

timestamps = sample(:, 1);
voltages = sample(:, 2);

Fsr = round(1/mean(diff(timestamps)))

Fsr = 200000

[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

ans = 0.1649

rms(voltages)   % rms of signal

ans = 0.3376

% 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]

Thank you very much. This is really helpful!! Please let me know how should I scale the ampltiude to a dB/V scaling factor of about 230 for this plot?

% %% PSD in terms of NON-DIMENSIONALIZed frequency
St= ((freq*Dj*0.001)/(Uj));
ss=(psd*Dj/Uj)/St;
plot(St, 10*log(ss))
title('Sound Pressure Level (SPL)');
xlabel('St');
ylabel('SPL (dB)');
xlim([0.1 1.5])

[Star Strider]

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.

.