Power spectral density plot converted to sound pressure level

Question

S_saleem 2023-6-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1977679-power-spectral-density-plot-converted-to-sound-pressure-level

评论： Star Strider 2023-6-4

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Hiro Yoshino 2023-6-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1977679-power-spectral-density-plot-converted-to-sound-pressure-level#answer_1249599

在 MATLAB Online 中打开

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

S_saleem 2023-6-3

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 2023-6-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1977679-power-spectral-density-plot-converted-to-sound-pressure-level#answer_1249614

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.

.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Star Strider 2023-6-3

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.

.

Star Strider 2023-6-3

编辑：Star Strider 2023-6-3

在 MATLAB Online 中打开

data..zip

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 2023-6-4

在 MATLAB Online 中打开

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 2023-6-4

在 MATLAB Online 中打开

data..zip

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 2023-6-4

在 MATLAB Online 中打开

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 2023-6-4

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.

.

请先登录，再进行评论。

Power spectral density plot converted to sound pressure level

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Power spectral density plot converted to sound pressure level

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

8 个评论 显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论