How to determine sampling frequency of wgn?

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

[Paul]

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]

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]

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

[William Rose]

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

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

[Paul]

@William Rose

I am confused because I thought the factor or 2 in the numerator for E(Pss(k)) only applies when computing the so-called 'one-sided' PSD, which is the default for periodogram. But the plot in the answer shows the frequency axis running to 10 Hz, which would be the two-sided periodogram.

So I reran the code from the answer and, as expected, the frequency axis runs to 5 Hz.

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

ans = 5

Pxxmn=mean(Pxx);
plot(f,Pxxmn,'-r.'); 
grid on; xlabel('Frequency (Hz)'); ylabel('Power/Hz'); title('P.S.D.')

[William Rose]

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

@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.