-
Notifications
You must be signed in to change notification settings - Fork 25
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Check white noise generation #89
Comments
The current formula to 'normalize' or 'scale' the standard deviation of the Gaussian distribution used to generate noise seems to have inconsistent units. Also, the from the derivation in the paper I gather that there was a disagreement/misunderstanding about what PSD plot of Gaussian white noise (GWN) with mean=0 and var=4. The curve represents the average over 2048 trials of 4 seconds. The Nyquist frequency (fmax) is 512Hz (fs = 1/dt and fmax = 1/2 fs). The PSD plot of Gaussian white noise with mean=0 and var=4. The curve represents the average over 2048 trials of 4 seconds. The Nyquist frequency (fmax) is 128Hz (fs = 1/dt and fmax = 1/2 fs). The shaded area represents the power between 0-100Hz. This value is constant regardless the value of Now, if we wish to keep the total power constant when the bandwidth changes (the shaded area would cover from -fmax to fmax), then the PSD (sigma**2) has to be adjusted and it will decrease proportionally to |
How do you calculate the PSD/Hz? I get a different result in Matlab with P128 = [];
P512 = [];
for j = 2048:-1:1
[P128(:,j),f128] = pwelch(2*randn(4*128,1),[],[],[],128);
[P512(:,j),f512] = pwelch(2*randn(4*512,1),[],[],[],512);
end
plot(f128,mean(P128,2))
hold on
plot(f512,mean(P512,2))
xlabel('Frequency')
ylabel('Power/Hz') Here it's the total power (integrated over all frequencies) that is the same |
Thank you for the example. Indeed, the two plots above indicate that the power Also, sometimes it's difficult to know whether people are talking about power If the variance of the Gaussian distribution is the same in the two cases, Is Parseval's theorem satisfied in the example above? Here is what I did: rng(42)
mu = 0; %mean of each gaussian process
L = 2048; %trials
T = 4; % s
dt = 2^-8; % s
t = dt:dt:T; % s
N = length(t); % samples
dt = t(2) - t(1); % s
fs = 1 / dt; % Hz
sigma = 2; % standard deviation (au/sqrt(Hz))
df = 1 / T; % Hz
if L == 1
z = mu + sigma*randn(1,N);
else
MU = mu*ones(1,N); % mean for all realizations
Cxx = (sigma^2)*diag(ones(N,1)); % cov mat
R = chol(Cxx);
% signals
z = repmat(MU,L,1) + randn(L,N)*R;
end
%%
Z = fft(z,[],2) / sqrt(N); %Scaling by sqrt(N);
if size(Z, 1) == 1
Pzavg = Z.*conj(Z);
else
Pzavg = mean(Z.*conj(Z)); %mean power spectral density
end
% Check Parseval's theorem
% frequency
P = norm(Z(1, :)).^2;
% time
p = norm(z(1, :)).^2;
%% Rearrange
f = (-N/2:N/2-1) * df;
Pzavg = fftshift(Pzavg); %Shift zero-frequency component to center of spectrum
%% Plot lines
subplot(1,1,1)
plot(f,Pzavg,'r');
axis([-fs/2 fs/2 0 10]); grid on;
ylabel('PSD(a.u.^{2}/Hz)');
xlabel('f[Hz]');
title({['PSD of Gaussian white noise \mu=0, \sigma^2=' num2str(sigma^2) ','],['f_{max}=' num2str(fs/2) '[Hz], \Delta f=' num2str(df) '[Hz],'], ['\Delta t=' num2str(dt) '[s] , T=' num2str(T) '[s]']}); I'm aware that the way power spectral density is computed for a given signal is In particular, The main issue is that it is difficult to properly justify the scaling of sigma by dt, dx, dy as it is in the code and to write a step-wise derivation of the scaling expression. |
Yes, Parseval's theorem is satisfied, at least as far as the definitions in Chris's notes. I just used
with the I'm admittedly not so comfortable in the continuous case - but there are at least a couple of indications that the variance of the process is expected to blow up as you increase the bandwidth. That's the result derived here. There's also a bit of discussion about the notation (second response, by Matt L). And see also the second post here. The other thing I'm not sure about is where you're using the norm to compute the energy - from here they associate the L2 norm with the energy of the signal, but both of those expressions are integrals so there's a The other question is about the |
Thank you for discussing this with me and sorry for the delay in replying (I was sick last week). And yes, I read Chris's notes and I see that it fits with your results and the examples are clear. About using the L2 norm: the square of the L2 norm gives the power of the signal and it should be the same of in the time domain and frequency domain. About _v(t)^2 = c: This is something that I didn't have very clear because I've read conflicting material. For instance, here they say: "The average value of the square of the signal (which is equal to the variance if the signal has zero mean) is equal to the integral of the power spectral density. " Ok. From that I understand the variance is the total power (PSD integrated over frequencies). Now, if you look at the first row of figure 3.1 you can see the timeseries, the power spectrum (psd I'm guessing) and autocorrelation. From these plots, I estimate that the variance used in the example is 1 (look at the autocorrelation at lag 0). But when I look at the height of the power spectrum which is 1 I have trouble believing that the power of the signal (PSD integrated over all frequencies) is equal to the variance. It's clearly not. I'm not sure what to make out of that. The integral over a fixed period of time is then the same irrespective of the sampling rate i.e. the energy depends on the variance and the duration of the signal, but not the sampling rate. An increase in the signal sampling rate effectively increases the bandwidth. In the case of white noise, because it has a constant power spectral density it means that now you have power from additional frequency components. That reasoning is only valid assuming that The factor of 1/sqrt(N) is used to have symmetry between the forward and inverse Fourier transforms |
No worries, hope you're feeling better! :) I think we're closing in on it. IMO the critical question boils down to What should change if you increase the bandwidth is the total power (integrated over all frequencies), while the variance (or alternatively PSD) remains constant. That's the part I can't reconcile with the discrete FFT normalisation - I think we agree that as the bandwidth increases, if the PSD (au^2/Hz) is kept the same, then the total power (integrated over frequency) is increased. But my understanding is that in the time domain, this corresponds to the variance increasing. What your calculation is saying is that if you sample the white noise at a higher sampling rate, the energy in the time domain increases. This is then related to the computation of the L2 norm: the square of the L2 norm gives the power of the signal and it should be the same of in the time domain and frequency domain. This is another crucial part - my understanding (from sources like these) is that this is an integral. The other reason I think this is the case - suppose we have a deterministic signal like a sine, and suppose it has a duration of 1s and a frequency of 1Hz (so we are sampling 1 cycle). If we sample it at 100Hz, there'll be 100 data points, and you could compute the power based on the signal norm. Now suppose it's sampled at 200Hz - now there are 200 data points, and we could compute the power based on the signal norm again. Shouldn't the energy in this signal should be the same in both cases? This is only true if the norm is computed with an integral. In the time domain, if the L2 norm is computed with an integral, then for white noise the power/energy won't go up if the sampling rate is increased. But without the integral, the vector norm will increase - I suspect this is the origin of our different results - Wolfram draws a distinction between the l2-norm (summation) and the L2-norm (integration), and I think it's the latter that refers to signal energy. I think once we agree on this, then we'll automatically agree on everything else! An increase in the signal sampling rate effectively increases the bandwidth. In the case of white noise, because it has a constant power spectral density it means that now you have power from additional frequency components. So I agree with this, but I would say, if we start in the frequency domain and increase the bandwidth, what happens to the autocorrelation? The autocorrelation is only nonzero at zero lag, so the only thing that can change is value of the autocorrelation at tau=0. So if we increase the bandwidth with constant PSD, then I think the autocorrelation at tau=0 increases. And then I would say that the autocorrelation function for a Gaussian process corresponds directly to the variance, so this implies the variance has increased. I think where we differ then is that you would say the autocorrelation of the Gaussian process requires an additional term that accounts for the sampling rate...? Another way I think about it is to consider synthesising a spectrally white time series from Fourier components, by setting all the components to the same amplitude and randomising their phases. Something like Then at a single point in time, the signal value is given by adding together N sines, each with a random phase. What are the statistics of this summation i.e. the statistics of the signal in the time domain? The mean should be zero, and the variance will depend on A_n (the PSD - the same for all components) and will increase with N. Increasing the bandwidth corresponds to having more terms, so this is also saying that as the bandwidth increases, so does the variance. Do you know what I'm missing here? Regarding that Figure 3.1, I don't trust my ability to interpret that figure, because it's not clear whether they are plotting power spectral density or not. Also, without units on their time axis, we can't tell how their bandwidth for the power spectrum came about, so it's possible that the spectrum was truncated. Their frequency spectrum x-axis doesn't have units either. Maybe it's just meant to be a qualitative figure |
Hey Guys,
Just a couple of comments that might be helpful.
I think that you are definitely correct that the FFT should be correctly scaled so that it approximates the integral.
However there seem to be some confusions arising out of Chris’s notes, at least taking them on first appearances. In his section 5, he states that the standard deviation of the white noise is sigma^2. By Parseval’s Theorem, sigma^2 = Power Spectrum Density x bandwidth.
If you want to implement Brownian Motion you want to take the bandwidth really large, and the power spectrum density constant. In other words, the standard deviation of your white noise should asymptote to infinity.
It seemed that in your notes, you are trying to keep the standard deviation to be constant, and asymptote the bandwidth to infinity and the power spectral density to zero. If you do this, then you should obtain a deterministic system in the limit of infinite bandwidth – i.e. as if there is no noise at all! I doubt in fact that this is what you really want from your model. Furthermore you are actually changing the physical model as you vary the bandwidth.
Good luck =)
James
From: Romesh Abeysuriya <[email protected]<mailto:[email protected]>>
Reply-To: BrainDynamicsUSYD/neurofield <[email protected]<mailto:[email protected]>>
Date: Tuesday, 30 May 2017 8:22 pm
To: BrainDynamicsUSYD/neurofield <[email protected]<mailto:[email protected]>>
Cc: Subscribed <[email protected]<mailto:[email protected]>>
Subject: Re: [BrainDynamicsUSYD/neurofield] Check white noise generation (#89)
No worries, hope you're feeling better! :) I think we're closing in on it. IMO the critical question boils down to
What should change if you increase the bandwidth is the total power (integrated over all frequencies), while the variance (or alternatively PSD) remains constant.
That's the part I can't reconcile with the discrete FFT normalisation - I think we agree that as the bandwidth increases, if the PSD (au^2/Hz) is kept the same, then the total power (integrated over frequency) is increased. But my understanding is that in the time domain, this corresponds to the variance increasing. What your calculation is saying is that if you sample the white noise at a higher sampling rate, the energy in the time domain increases. This is then related to the computation of the L2 norm:
the square of the L2 norm gives the power of the signal and it should be the same of in the time domain and frequency domain.
This is another crucial part - my understanding (from sources<http://users.abo.fi/htoivone/courses/robust/rob2.pdf> like<https://math.stackexchange.com/questions/885998/physical-interpretation-of-l-1-and-l-2-norms> these<http://uqlab.github.io/pdfs/aero632/02%20Signal%20Norms.pdf>) is that this is an integral. The other reason I think this is the case - suppose we have a deterministic signal like a sine, and suppose it has a duration of 1s and a frequency of 1Hz (so we are sampling 1 cycle). If we sample it at 100Hz, there'll be 100 data points, and you could compute the power based on the signal norm. Now suppose it's sampled at 200Hz - now there are 200 data points, and we could compute the power based on the signal norm again. Shouldn't the energy in this signal should be the same in both cases? This is only true if the norm is computed with an integral.
In the time domain, if the L2 norm is computed with an integral, then for white noise the power/energy won't go up if the sampling rate is increased. But without the integral, the vector norm will increase - I suspect this is the origin of our different results - Wolfram<http://mathworld.wolfram.com/L2-Norm.html> draws a distinction between the l2-norm (summation) and the L2-norm (integration), and I think it's the latter that refers to signal energy.
I think once we agree on this, then we'll automatically agree on everything else!
An increase in the signal sampling rate effectively increases the bandwidth. In the case of white noise, because it has a constant power spectral density it means that now you have power from additional frequency components.
So I agree with this, but I would say, if we start in the frequency domain and increase the bandwidth, what happens to the autocorrelation? The autocorrelation is only nonzero at zero lag, so the only thing that can change is value of the autocorrelation at tau=0. So if we increase the bandwidth with constant PSD, then I think the autocorrelation at tau=0 increases. And then I would say that the autocorrelation function for a Gaussian process corresponds directly to the variance, so this implies the variance has increased. I think where we differ then is that you would say the autocorrelation of the Gaussian process requires an additional term that accounts for the sampling rate...?
Another way I think about it is to consider synthesising a spectrally white time series from Fourier components, by setting all the components to the same amplitude and randomising their phases. Something like
[latex-image-3]<https://cloud.githubusercontent.com/assets/755796/26578934/b0318852-4529-11e7-8252-c9b45af3f7f7.png>
Then at a single point in time, the signal value is given by adding together N sines, each with a random phase. What are the statistics of this summation i.e. the statistics of the signal in the time domain? The mean should be zero, and the variance will depend on A_n (the PSD - the same for all components) and will increase with N. Increasing the bandwidth corresponds to having more terms, so this is also saying that as the bandwidth increases, so does the variance. Do you know what I'm missing here?
Regarding that Figure 3.1, I don't trust my ability to interpret that figure, because it's not clear whether they are plotting power spectral density or not. Also, without units on their time axis, we can't tell how their bandwidth for the power spectrum came about, so it's possible that the spectrum was truncated. Their frequency spectrum x-axis doesn't have units either. Maybe it's just meant to be a qualitative figure
—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub<#89 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AM9RjpSekprCP0Wz-r1QUsZBbsVXFMUeks5r--30gaJpZM4MFy1L>.
|
Thanks James, I was hoping you'd join in :) It seemed that in your notes, you are trying to keep the standard deviation to be constant, and asymptote the bandwidth to infinity and the power spectral density to zero. If you do this, then you should obtain a deterministic system in the limit of infinite bandwidth – i.e. as if there is no noise at all! I doubt in fact that this is what you really want from your model. Yes, that's indeed not what I meant. The issue is that for the linear analytic spectrum, the PSD at say, 10Hz, is fixed, and of course you can evaluate the spectrum/transfer function at whatever set of frequencies you like because they are all independent. But when we integrate numerically, how do we ensure that the PSD of the noise in the bin spanning 10Hz has the same value that went into the theoretical calculation? So the aim is to keep the PSD constant as the bandwidth (in both frequency and wavenumber) changes, by adjusting the variance accordingly. The current scaling in NeuroField does what you've just said, where as the time step approaches 0 and the bandwidth thus approaches infinity, the variance also approaches infinity. So at the very least, the clarity of the derivation needs to be improved! |
Lol yes well I must seem very opinionated to you guys ;) Just can’t help commenting sometimes…
Yes, that's indeed not what I meant. The issue is that for the linear analytic spectrum, the PSD at say, 10Hz, is fixed, and of course you can evaluate the spectrum/transfer function at whatever set of frequencies you like because they are all independent. But when we integrate numerically, how do we ensure that the PSD of the noise in the bin spanning 10Hz has the same value that went into the theoretical calculation? So the aim is to keep the PSD constant as the bandwidth (in both frequency and wavenumber) changes, by adjusting the variance accordingly. The current scaling in NeuroField does what you've just said, where as the time step approaches 0 and the bandwidth thus approaches infinity, the variance also approaches infinity.
I agree with what you have written above in principle. I presume that you scale the variance to be proportional to the bandwidth of the power spectrum. It might be that I misunderstood what was written in the notes, or that the version that I read was outdated.
As a suggestion, a quick way of checking that you have scaled the white noise correctly is to check that the order of magnitude of your solution is approximately invariant under a change of the timestep / bandwidth. That is, take the timestep to be really small, run the simulation, and then halve the timestep and run it again. Obviously the solutions will be different because it is stochastic, but the orders of magnitude of the solutions should be approximately the same. On first appearances, in the formula I saw for how you calculate the variance of the white noise, it seemed that the order of magnitude of the solution will change as you change the scaling of the timestep, which would be a red flag that something is wrong.
In any case I am sure that you don’t need me to tell you what to do …. good luck =)
|
I don't think anything has changed since this document, is that the one you're talking about? There, the equation at the top of page 4 has the variance being linearly proportional to the bandwidth (via the sampling rate) when the PSD (\phi^2_n(\omega)) is held constant [edit: note that this is the PSD that goes into the analytic calculation, which is selected independently in advance]. How should the magnitude of the solution be measured? In the time domain or in the frequency domain? I would expect that the power in a specific frequency range e.g. 1-100Hz should remain the same, but if the PSD has a fixed value and the bandwidth is increased, then the total power in the frequency domain is higher, and thus the power in the time domain should be higher as well i.e. the magnitude of the solution should grow with the variance, in the absence of any filtering, which reflects the definition of the white noise...? |
I have just glanced at this. There are things that I do not understand on first appearances. Maybe there are some typos?
On the top of page 3 you seem to write that the autocorrelation function for white noise is equal to sigma^2
But you then write that the power spectral density is sigma^2? In your previous email, you agreed with me that the variance of white noise is equal to power spectral noise times bandwidth….
Which would imply that sigma^2 = sigma^2 x bandwidth? This is a contradiction unless bandwidth=1.
Are we using the same terminology? The variance is by definition the value of the autocorrelation function at tau=0, i.e. the expectation (average over all realizations) of the white noise signal. It should be very very large for a white noise process.
From: Romesh Abeysuriya <[email protected]<mailto:[email protected]>>
Reply-To: BrainDynamicsUSYD/neurofield <[email protected]<mailto:[email protected]>>
Date: Friday, 2 June 2017 6:03 pm
To: BrainDynamicsUSYD/neurofield <[email protected]<mailto:[email protected]>>
Cc: James MacLaurin <[email protected]<mailto:[email protected]>>, Comment <[email protected]<mailto:[email protected]>>
Subject: Re: [BrainDynamicsUSYD/neurofield] Check white noise generation (#89)
I don't think anything has changed since this document<https://github.com/BrainDynamicsUSYD/neurofield/blob/master/doc/noise.pdf>, is that the one you're talking about? There, the equation at the top of page 4 has the variance being linearly proportional to the bandwidth (via the sampling rate) when the PSD (\phi^2_n(\omega)) is held constant.
How should the magnitude of the solution be measured? In the time domain or in the frequency domain? I would expect that the power in a specific frequency range e.g. 1-100Hz should remain the same, but if the PSD has a fixed value and the bandwidth is increased, then the total power in the frequency domain is higher, and thus the power in the time domain should be higher as well i.e. the magnitude of the solution should grow with the variance, in the absence of any filtering, which reflects the definition of the white noise...?
—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub<#89 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AM9Rjlb8IvUn2tDKja83WlquVVwYp9wgks5r_8HngaJpZM4MFy1L>.
|
Yes I think the top of Page 3 is incorrect because there I was quoting the textbook formulae that have the issues with notation discussed in the second response here. I knew there was a problem, but at the time I attributed this to the result corresponding to an ideal continuous case, rather than being an issue with the notation. So those formulae aren't actually used in the derivation of the result that's in NeuroField :) The scaling originates from the third equation on that page, where the expectation value of the white noise process (i.e. its variance) is equated to the total power in the spectrum integrated over frequency (i.e. accounting for bandwidth) |
Hi Romesh,
I mostly agree with the scalings in the rest of your notes. Sometimes there are factors of pi that I wouldn’t really expect, but aside from this it looks ok.
Good luck with the paper =)
James
|
The current implementation of White noise assumes we're working with a square grid in order to compute
deltax
. If the grid is rectangular, the current scaling of the noise will result in an incorrect value of deltax.https://github.com/BrainDynamicsUSYD/nftsim/blob/master/src/timeseries.cpp#L330
Related issues #90 and #151.
Current behaviour if the grid has 450 nodes and Longside is 30 (nodes), then
Changing the time step to 4.8828e-04 and repeating the analysis we obtain:
Bottom line: the scaling of the variance of the underlying normal distribution used to generate a sequence of random numbers, really depends on what property of the resulting signal one needs to keep constant (eg, PSD or average power). This scaling might be at odds with interpretations of what that sequence of random numbers is (eg, increments of a Wiener process, just and input signal with a really broad band), etc.
The text was updated successfully, but these errors were encountered: