How to determine sampling frequency of wgn?

Question

Jan 2024-7-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2134261-how-to-determine-sampling-frequency-of-wgn

评论： William Rose 2024-7-5

采纳的回答： William Rose

在 MATLAB Online 中打开

Hello there,

I'm trying to understand noise analysis and the concept of power spectral density. I've found this article:

https://pure.mpg.de/pubman/faces/ViewItemOverviewPage.jsp?itemId=item_152164

Which was really helpful and interesting, but I've encountered some problems while trying to put it in practise.

For convenience I wanted to generate white noise (mainly because of its flat spectral density) via wgn funtion and then calculate its PSD. If i understood it correctly, wgn creates discrete samples, which completely lack any relatationship to time.

And here comes my problem - while normalizing the results, equivalent noise bandwidth is dependent on sampling rate and so is PSD:

,

where w are window function values, S are complex results from fft function.

By experimenting with code I came upon sampling rate of 2 samples/s, when mean power of PSD was getting close to 5 dB. Unfortunatelly, I still cannot confidently explain this.

Is there please any way to obtain "realistic" white noise samples in MATLAB? Or do I have to add white noise to (co)sine signal with known sampling? If I'm missing something, please let me know.

Thank you in advance,

Jan

CODE:

% Size parameters

L = pow2(16);

iterations = 100;

% ??? (Guessed value)

fs = 2;

% Frequency resolution

f_res = fs/L;

% Create noise with power of 5 dBW

data = wgn(L,iterations, 5);

% Apply Hann window

h = hann(L);

w = ones(L,iterations) .* h;

data = data .* w;

% Calculate normalizing factor S2

S2 = sum(h.^2);

% Calculate fft

S = fft(data,L,1);

S = S .* conj(S);

% Average power spectrum

for i = 1:L

PSD(i) = mean(S(i,:));

end

% Normalize results

PSD = 2 * PSD / (fs * S2);

% Discard symmectric part

PSD = PSD(1:(L/2+1));

% PSD in dB(W)

PSD = 10*log10(PSD);

% Mean value of PSD

display("Power = " + mean(PSD) + " dB");

"Power = 4.9725 dB"

% Calculate frequency

f = (0:(L/2)) * f_res;

% Plot data

figure

plot(f,PSD)

xlabel("f [Hz]")

ylabel("PSD [dB/Hz]")

2 个评论
显示无隐藏无

Paul 2024-7-4

编辑：Paul 2024-7-4

Hi Jan,

Can you clarify what the question? What exactly are you not confident in explaining and what is mean by "realistic" white noise?

Jan 2024-7-4

在 MATLAB Online 中打开

Hi @Paul,

I apologize for the confusion I've made. When measuring noise in real-life scenarion, each iteration would result in two vectors - first containing signal magnitude and the second one time stamps. Or just a magnitude vector and value of sampling rate. Either way, I would be able to normalize PSD correctly, since I could calculate/would know the

.

So, my question was, if there is a way to simulate such a measurement in MATLAB, that would keep its relationship to time or sampling? Something like:

[mag, time] = wgn(L, iterations, power, fs);

but as I said before, I know wgn creates discrete samples. I just need to find a way to create white noise dependent on sampling rate, so when

changes, the mean value of normalized PSD reamains at the approximately same level.

请先登录，再进行评论。

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

Answer 1

William Rose 2024-7-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2134261-how-to-determine-sampling-frequency-of-wgn#answer_1481296

编辑：William Rose 2024-7-5

在 MATLAB Online 中打开

@Jan,

[Edit: Re-ran the script because the plot below was wrong, probably from an earlier version of the script. Thanks to @Paul for catching it. Plot looks right now.]

You wrote "I just need to find a way to create white noise dependent on sampling rate, so when changes, the mean value of normalized PSD reamains at the approximately same level."

The expected power spectral density (PSD) associated with discretely sampled Gaussian white noise is

where fs=sampling rate, and Var(x)=variance of the time domain signal. Note that this is expected PSD. The actual PSD will have random fluctuation, To check if the formula above is correct, let's compute the PSD of an N point signal, M times, and average the M PSDs.

M=1000; N=256;

s=1; % standard deviation of the noise

fs=10; % sampling rate (Hz)

Pxx=zeros(M,N/2+1); % allocate array to hold the PSD results

for i=1:M

x=s*randn(1,N);

[Pxx(i,:),f]=periodogram(x,[],N,fs);

end

Pxxmn=mean(Pxx);

plot(f,Pxxmn,'-r.');

grid on; xlabel('Frequency (Hz)'); ylabel('Power/Hz'); title('P.S.D.')

The equation above predicts that the expected PSD should be 2*(1^2)/10=0.2. This is what we observe in the plot above.

Therefore, if you want to add noise so that its PSD is independent of the sampling rate, you need to scale the variance of the noise in proportion to the sampling rate. For example, if you increase fs by a factor of 2, then you should increase the standard deviation of the added noise by sqrt(2).

5 个评论
显示 3更早的评论隐藏 3更早的评论

William Rose 2024-7-4

编辑：William Rose 2024-7-4

在 MATLAB Online 中打开

@Jan,

[Edit: Change "the PSDs ... are the same" to "the mean PSDs ... are the same (except for the expected random fluctuations)", in my remarks.]

Example: Create two random signals, each of duration T=2 seconds. Let fs=50 Hz for signal x1 and fs=100 Hz for x2. The desired expected PSD (power per Hz) for both signals is 10. Plot the mean PSD (averaged over M=1000 interations) for x1 and x2.

The script below uses the formula

where s=standard deviation of x(t), fs=sampling rate, and PSDdes=desired PSD. This formula follows from the equation in my answer above.

% Define constants

fs1=50; fs2=100; % sampling rates (Hz)

T=2; % signal duration (s)

M=1000; % number of iterations

PSDdes=10; % desired PSD (power/Hz)

% Compute signals and PSDs

s1=sqrt(PSDdes*fs1/2); % std.dev. of x1

s2=sqrt(PSDdes*fs2/2); % std.dev. of x2

N1=T*fs1; N2=T*fs2; % signal duration (points)

Pxx1=zeros(M,N1/2+1); % allocate array to hold PSDs for x1

Pxx2=zeros(M,N2/2+1); % allocate array to hold PSDs for x2

for i=1:M

x1=s1*randn(1,N1);

[Pxx1(i,:),f1]=periodogram(x1,[],N1,fs1);

x2=s2*randn(1,N2);

[Pxx2(i,:),f2]=periodogram(x2,[],N2,fs2);

end

Pxx1mn=mean(Pxx1);

Pxx2mn=mean(Pxx2);

% Plot results

plot(f1,Pxx1mn,'-r.',f2,Pxx2mn,'-b.');

grid on; xlabel('Frequency (Hz)'); ylabel('Power/Hz')

legend('x1','x2'); title('P.S.D.')

The mean PSDs for x1 and x2 are the same (except for the expected random fluctuations), and have the desired value (=10) in the frequency range where the PSD exists for each signal. Of course the PSD goes up to 25 Hz for x1 (with fs1=50 Hz), and the PSD goes up to 50 Hz for x2 (with fs2=100 Hz).

Jan 2024-7-4

在 MATLAB Online 中打开

Dear @William Rose,

thank you very much. I didn't know about such a relation for expected power. I believe it solves all my problems, including why I could find expected power of 5 dB only at 2 Hz sampling rate. I've edited my code, thinking that maybe someone will find it useful:

% Data parameters

L = pow2(16);

iterations = 100;

% Sampling rate

fs = 10000;

% Frequency resolution

f_res = fs/L;

% Create noise

data = wgn(L,iterations, 5) * sqrt(fs/2);

% Estimate power

for i = 1:iterations

P_it(i) = 2 * var(data(:,i)) / fs;

end

P = 10*log10(mean(P_it));

display("Power estimate: " + P + " dB");

"Power estimate: 4.9992 dB"

% Apply Hann window

h = hann(L);

win = ones(L,iterations) .* h;

data = data .* win;

% Calculate normalizing factor S2

S2 = sum(h.^2);

% Calculate fft

S = fft(data,L,1);

S = S .* conj(S);

% Average power spectrum

for i = 1:L

PSD(i) = mean(S(i,:));

end

% Normalize results

PSD = 2 * PSD / (fs * S2);

% Discard symmectric part

PSD = PSD(1:(L/2+1));

% PSD in dB

PSD = 10*log10(PSD);

display("Average PSD: " + mean(PSD) + " dB");

"Average PSD: 4.9765 dB"

% Calculate frequency

f = (0:(L/2)) * f_res;

% Plot data

figure

plot(f,PSD)

hold on

plot(f,ones(1,L/2+1)*P)

xlabel("f [Hz]")

ylabel("PSD [dB/Hz]")

legend("PSD", "Power Estimate")

Best wisches,

Jan

William Rose 2024-7-5

@Paul, Thank you for your attention to detail. Of course you are correct that the frequency axis should stop at fs/2, since periodogram returns the one-sided spectrum by default.

I suspect that I ran the script with fs=20 Hz, then changed fs to 10 Hz in the script, but forgot to re-run the script.

The plot with two PSDs, in my comment, shows the one-sided spectra extending up to frequency=fs/2, as expected.

William Rose 2024-7-5

@Jan,

You're welcome. You wrote "I didn't know about such a relation for expected power."

To be honest, I did not know anbout that relationship either, at least not in that exact form, until I worked it out from Parseval's theorem.

You wrote "I believe it solves all my problems..." I like your sense of humor.

请先登录，再进行评论。

How to determine sampling frequency of wgn?

2 个评论
显示无隐藏无

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How to determine sampling frequency of wgn?

2 个评论 显示 无隐藏 无

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

5 个评论
显示 3更早的评论隐藏 3更早的评论