Prev: Image rotation detection probelm
Next: Filling in the space between two parallel lines with a color
From: Dean on 2 Aug 2010 21:07
Hi guys,

I understand that a PSD of a white Gaussian noise of is equivalent to its variance. However, I can't even close to get this value using Periodogram routine. I believe I've forget a scale factor somewhere.

For example, I've a WGN of variance of 16, I run a very large number of data sequence through Periodogram, and do a linear fit to the spectrum. Although the PSD looks flat but the value is wrong.

According to Parseval's Theorem, the area under the spectrum must equal the variance and therefore equivalent to PSD. I found using numerical integration of the area gives me a better estimation.

On the other hand, I expect the variance is also equivalent to the square of the power spectrum. In which I failed to get a good estimation.

Anyone knows how to correctly estimate the variance or has any comments on this please let me know.

Cheers,

Dean
--------------------Code Here------------------------------------------------
N = 100000;
H = spectrum.periodogram;
x = 4*randn(N,1); % wgn with a variance of 16
Py = psd(H,x);
% in dB PSD is flat and equal to 10*log10(16)
plot(Py.Frequencies,10*log10(Py.Data));
% then do a linear fit to check if the PSD_hat=10*log10(16), sadly not as I expected.
trapz(Py.Data)/N % this seem works
sumsqr(Py.Data)/N % this one doesn't
From: Greg Heath on 3 Aug 2010 05:42
On Aug 2, 9:07 pm, "Dean " <jiangwei0...(a)gmail.com> wrote:
> Hi guys,
>
> I understand that a PSD of a white Gaussian noise of is equivalent to its variance. However, I can't even close to get this value using Periodogram routine. I believe I've forget a scale factor somewhere.
>
> For example, I've a WGN of variance of 16, I run a very large number of data sequence through Periodogram, and do a linear fit to the spectrum. Although the PSD looks flat but the value is wrong.
>
> According to Parseval's Theorem, the area under the spectrum must equal the variance and therefore equivalent to PSD. I found using numerical integration of the area gives me a better estimation.
>
> On the other hand, I expect the variance is also equivalent to the square of the power spectrum. In which I failed to get a good estimation.
>
> Anyone knows how to correctly estimate the variance or has any comments on this please let me know.
>
> Cheers,
>
> Dean
> --------------------Code Here------------------------------------------------
> N = 100000;
> H = spectrum.periodogram;
> x = 4*randn(N,1); % wgn with a variance of 16
> Py = psd(H,x);
> % in dB PSD is flat and equal to 10*log10(16)
> plot(Py.Frequencies,10*log10(Py.Data));
> % then do a linear fit to check if the PSD_hat=10*log10(16), sadly not as I expected.
> trapz(Py.Data)/N % this seem works
> sumsqr(Py.Data)/N % this one doesn't

close all, clear all, clc
N = 100000; % N = 100 might be better
randn('state',0)
x = 4*randn(N,1); % wgn with a variance of 16
meanx = mean(x) % 0.012829
varx = var(x) % 16.118
meansqrx = mean(x.^2) % 16.118

% NOTE: meanx is to small and N is too large to
% distinguish between the two

check = meansqrx-varx % 3.4184e-006

X = fft(x);
PSD = abs(X).^2/N; % Not scaled by (1/Fs)
meanPSD = mean(PSD) % 16.118
PSDdb = 10*log10(PSD);
meanPSDdb = mean(PSDdb) % 9.5829
illegal = 10*log10(meanPSD) % 12.073

Hope this helps.

Greg
From: Dean on 4 Aug 2010 04:04
Greg Heath <heath(a)alumni.brown.edu> wrote in message <2ba31fe8-f154-4dec-ad32-64cbae041191(a)w12g2000yqj.googlegroups.com>...
> On Aug 2, 9:07 pm, "Dean " <jiangwei0...(a)gmail.com> wrote:
> > Hi guys,
> >
> > I understand that a PSD of a white Gaussian noise of is equivalent to its variance. However, I can't even close to get this value using Periodogram routine. I believe I've forget a scale factor somewhere.
> >
> > For example, I've a WGN of variance of 16, I run a very large number of data sequence through Periodogram, and do a linear fit to the spectrum. Although the PSD looks flat but the value is wrong.
> >
> > According to Parseval's Theorem, the area under the spectrum must equal the variance and therefore equivalent to PSD. I found using numerical integration of the area gives me a better estimation.
> >
> > On the other hand, I expect the variance is also equivalent to the square of the power spectrum. In which I failed to get a good estimation.
> >
> > Anyone knows how to correctly estimate the variance or has any comments on this please let me know.
> >
> > Cheers,
> >
> > Dean
> > --------------------Code Here------------------------------------------------
> > N = 100000;
> > H = spectrum.periodogram;
> > x = 4*randn(N,1); % wgn with a variance of 16
> > Py = psd(H,x);
> > % in dB PSD is flat and equal to 10*log10(16)
> > plot(Py.Frequencies,10*log10(Py.Data));
> > % then do a linear fit to check if the PSD_hat=10*log10(16), sadly not as I expected.
> > trapz(Py.Data)/N % this seem works
> > sumsqr(Py.Data)/N % this one doesn't
>
> close all, clear all, clc
> N = 100000; % N = 100 might be better
> randn('state',0)
> x = 4*randn(N,1); % wgn with a variance of 16
> meanx = mean(x) % 0.012829
> varx = var(x) % 16.118
> meansqrx = mean(x.^2) % 16.118
>
> % NOTE: meanx is to small and N is too large to
> % distinguish between the two
>
> check = meansqrx-varx % 3.4184e-006
>
> X = fft(x);
> PSD = abs(X).^2/N; % Not scaled by (1/Fs)
> meanPSD = mean(PSD) % 16.118
> PSDdb = 10*log10(PSD);
> meanPSDdb = mean(PSDdb) % 9.5829
> illegal = 10*log10(meanPSD) % 12.073
>
> Hope this helps.
>
> Greg

Hi Greg,

My main purpose is to find a relationship between variance and Periodogram PSD estimations. The reason I use very large value N is because I believe the value will converge to its true value.

For example, If I use 100, varx=12.0688 where the theoretical value should be 16.

I started with WGN because I know its variance and mean of PSD as N->infinity. I don't know any other "Gold standard", I thought it might be wise to start with something I know. (highly appreciate if anyone as any input/comments on this topic).

For a WGN with variance = 1, its mean PSD and variance are both = 1. This seems doesn't agree with your comments.

I am sorry for my understanding, and I am grateful for your answer.

Dean
From: Wayne King on 4 Aug 2010 07:00
"Dean " <jiangwei0515(a)gmail.com> wrote in message <i3b6tk$1a1$1(a)fred.mathworks.com>...
> Greg Heath <heath(a)alumni.brown.edu> wrote in message <2ba31fe8-f154-4dec-ad32-64cbae041191(a)w12g2000yqj.googlegroups.com>...
> > On Aug 2, 9:07 pm, "Dean " <jiangwei0...(a)gmail.com> wrote:
> > > Hi guys,
> > >
> > > I understand that a PSD of a white Gaussian noise of is equivalent to its variance. However, I can't even close to get this value using Periodogram routine. I believe I've forget a scale factor somewhere.
> > >
> > > For example, I've a WGN of variance of 16, I run a very large number of data sequence through Periodogram, and do a linear fit to the spectrum. Although the PSD looks flat but the value is wrong.
> > >
> > > According to Parseval's Theorem, the area under the spectrum must equal the variance and therefore equivalent to PSD. I found using numerical integration of the area gives me a better estimation.
> > >
> > > On the other hand, I expect the variance is also equivalent to the square of the power spectrum. In which I failed to get a good estimation.
> > >
> > > Anyone knows how to correctly estimate the variance or has any comments on this please let me know.
> > >
> > > Cheers,
> > >
> > > Dean
> > > --------------------Code Here------------------------------------------------
> > > N = 100000;
> > > H = spectrum.periodogram;
> > > x = 4*randn(N,1); % wgn with a variance of 16
> > > Py = psd(H,x);
> > > % in dB PSD is flat and equal to 10*log10(16)
> > > plot(Py.Frequencies,10*log10(Py.Data));
> > > % then do a linear fit to check if the PSD_hat=10*log10(16), sadly not as I expected.
> > > trapz(Py.Data)/N % this seem works
> > > sumsqr(Py.Data)/N % this one doesn't
> >
> > close all, clear all, clc
> > N = 100000; % N = 100 might be better
> > randn('state',0)
> > x = 4*randn(N,1); % wgn with a variance of 16
> > meanx = mean(x) % 0.012829
> > varx = var(x) % 16.118
> > meansqrx = mean(x.^2) % 16.118
> >
> > % NOTE: meanx is to small and N is too large to
> > % distinguish between the two
> >
> > check = meansqrx-varx % 3.4184e-006
> >
> > X = fft(x);
> > PSD = abs(X).^2/N; % Not scaled by (1/Fs)
> > meanPSD = mean(PSD) % 16.118
> > PSDdb = 10*log10(PSD);
> > meanPSDdb = mean(PSDdb) % 9.5829
> > illegal = 10*log10(meanPSD) % 12.073
> >
> > Hope this helps.
> >
> > Greg
>
> Hi Greg,
>
> My main purpose is to find a relationship between variance and Periodogram PSD estimations. The reason I use very large value N is because I believe the value will converge to its true value.
>
> For example, If I use 100, varx=12.0688 where the theoretical value should be 16.
>
> I started with WGN because I know its variance and mean of PSD as N->infinity. I don't know any other "Gold standard", I thought it might be wise to start with something I know. (highly appreciate if anyone as any input/comments on this topic).
>
> For a WGN with variance = 1, its mean PSD and variance are both = 1. This seems doesn't agree with your comments.
>
> I am sorry for my understanding, and I am grateful for your answer.
>
> Dean

Hi Dean, keep in mind that the unmodified periodogram only converges in expected value (the mean) to true PSD as N goes to infinity. The periodogram is no a consistent estimator of the true PSD, which means that its variance will not go to zero. The variance of the periodogram is approximately the square of the true PSD.

Wayne
From: Dean on 4 Aug 2010 07:26
> > Dean
>
> Hi Dean, keep in mind that the unmodified periodogram only converges in expected value (the mean) to true PSD as N goes to infinity. The periodogram is no a consistent estimator of the true PSD, which means that its variance will not go to zero. The variance of the periodogram is approximately the square of the true PSD.
>
> Wayne

Hi Wayne,

That's exactly I've been told by textbooks. I'm trying my best to verify it trough simulations but didn't succeed. Can you suggest a way to verify it?

Cheers,

Dean
 |  Next  |  Last
Pages: 1 2
Prev: Image rotation detection probelm
Next: Filling in the space between two parallel lines with a color