Convolution of same PSDs

Hi,

I have a (for simplicity brickwall) bandlimited white Gaussian noise
signal x[n]. Hence I know roughly its PSD.

I would like to get an estimate of ||x^3[n]||/||x[n]||^3  (||.|| l2-Norm
of the discrete time signal x[n]).

My idea was now to go to frequency domain (Parseval) and obtain the PSD
of x^3[n] via 2-fold convolution of X[k].

If I would just approximate the PSD of x[n] as ideal square function,
then X[k]*X[k]*X[k] would be a quadratic and I can find an analytic
estimate. Of course, this does not match at all with the numeric results
when I have a random signal.

Any hints on that?

Regards,
Peter

On Sat, 04 Apr 2015 01:01:59 -0700, Peter Mairhofer wrote:

> Hi,
> 
> I have a (for simplicity brickwall) bandlimited white Gaussian noise
> signal x[n]. Hence I know roughly its PSD.
> 
> I would like to get an estimate of ||x^3[n]||/||x[n]||^3  (||.|| l2-Norm
> of the discrete time signal x[n]).

Your notation isn't clear here.  Are you referring to x[n] as the n'th 
sample of x, here, or are you referring to x[n] as the whole signal?  If 
it's the "whole signal" and the signal is of infinite extent, I don't 
think the answer is finite.

> My idea was now to go to frequency domain (Parseval) and obtain the PSD
> of x^3[n] via 2-fold convolution of X[k].

This should be valid.

> If I would just approximate the PSD of x[n] as ideal square function,
> then X[k]*X[k]*X[k] would be a quadratic and I can find an analytic
> estimate. Of course, this does not match at all with the numeric results
> when I have a random signal.

Do you mean you've tried it and it doesn't work?  Are you sure you have 
enough samples?

If the expression you're trying to get estimates of means what I think it 
does, then you're dividing by a quantity with a nonzero probability of 
being zero, so your expression can go to infinity -- are you sure the 
math isn't just blowing up in your face?

-- 
www.wescottdesign.com

On 2015-04-04 7:57, Tim Wescott wrote:
> On Sat, 04 Apr 2015 01:01:59 -0700, Peter Mairhofer wrote:
> 
>> Hi,
>>
>> I have a (for simplicity brickwall) bandlimited white Gaussian noise
>> signal x[n]. Hence I know roughly its PSD.
>>
>> I would like to get an estimate of ||x^3[n]||/||x[n]||^3  (||.|| l2-Norm
>> of the discrete time signal x[n]).
> 
> Your notation isn't clear here.  Are you referring to x[n] as the n'th 
> sample of x, here, or are you referring to x[n] as the whole signal?  If 
> it's the "whole signal" and the signal is of infinite extent, I don't 
> think the answer is finite.

It is the whole signal.

For simplicity I assume now I have a signal of, say, length N=20400 with
the above properties.

One instance of the signal is generated with:

% WGN
x = sqrt(sigma2) * randn(N,1);

% BLWGN, brickwall filter
X = fft(x);
X(2041:18361) = 0;
x = ifft(X);

Mathematically I can just normalize by some period and let period go to
infinity in the limit.

Continuing with MATLAB syntax, I would seek an (analytical) estimate for
norm(x.^3)/norm(x)^3

>> My idea was now to go to frequency domain (Parseval) and obtain the PSD
>> of x^3[n] via 2-fold convolution of X[k].
> 
> This should be valid.

See below.

>> If I would just approximate the PSD of x[n] as ideal square function,
>> then X[k]*X[k]*X[k] would be a quadratic and I can find an analytic
>> estimate. Of course, this does not match at all with the numeric results
>> when I have a random signal.
> 
> Do you mean you've tried it and it doesn't work?  Are you sure you have 
> enough samples?

Not sure.

> If the expression you're trying to get estimates of means what I think it 
> does, then you're dividing by a quantity with a nonzero probability of 
> being zero, so your expression can go to infinity -- are you sure the 
> math isn't just blowing up in your face?

I don't think so. I don't get means but mean-squares, sum(x.^6) /
sum(x.^2)^3. Both expressions should not be near zero (it's essentially
"signal energy")

What I want is not to do the *numerical* estimate but derive an
*analytical* expression for the estimate based on the assumption that
x[n] is BLWGN described by zero mean, BW, variance.

I think my assumption is wrong to just replace my spectrum with the
analytical PSD. On the other hand I am very sure it should be possible
to derive an analytical expression for the estimate.

Let me quickly continue elaborating what I did:
Given the signal generated above, I always get a similar *numerical* answer:

>> 20*log10(norm(x.^3)/norm(x)^3)
ans =

  -74.3041

Of course, also if I do it with DFT (works due to Parseval):

>> X = fft(x);
>> X3 = fftshift(cconv(cconv(fftshift(X),fftshift(X),N),fftshift(X),N)/N^2);
>> 20*log10(norm(X3)/norm(X)^3*N)
ans =

  -74.3041

Now the "idea" from above was to replace "X" with the analytical form of
the signals PSD:

X = 0*X;
X(1:2040) = c;
X(18362:end) = c;
x = ifft(X);

But now (same result in DFT domain and convolution):

>> 20*log10(norm(x.^3)/norm(x)^3)
ans =

  -16.5779


So I can't just do that. Probably what you mean is I would need to
extend the signal in time somehow? Maybe there is a different approach
at all? If it's easier to approach it in the continuous-time domain I'm
also open for hints on that.


Thank you,
Peter

Let me add something below:

On 2015-04-04 12:34, Peter Mairhofer wrote:
> On 2015-04-04 7:57, Tim Wescott wrote:
>> On Sat, 04 Apr 2015 01:01:59 -0700, Peter Mairhofer wrote:
>>
>>> Hi,
>>>
>>> I have a (for simplicity brickwall) bandlimited white Gaussian noise
>>> signal x[n]. Hence I know roughly its PSD.
>>>
>>> I would like to get an estimate of ||x^3[n]||/||x[n]||^3  (||.|| l2-Norm
>>> of the discrete time signal x[n]).
>>
>> Your notation isn't clear here.  Are you referring to x[n] as the n'th 
>> sample of x, here, or are you referring to x[n] as the whole signal?  If 
>> it's the "whole signal" and the signal is of infinite extent, I don't 
>> think the answer is finite.
> 
> It is the whole signal.
> 
> For simplicity I assume now I have a signal of, say, length N=20400 with
> the above properties.
> 
> One instance of the signal is generated with:
> 
> % WGN
> x = sqrt(sigma2) * randn(N,1);
> 
> % BLWGN, brickwall filter
> X = fft(x);
> X(2041:18361) = 0;
> x = ifft(X);
> 
> Mathematically I can just normalize by some period and let period go to
> infinity in the limit.
> 
> Continuing with MATLAB syntax, I would seek an (analytical) estimate for
> norm(x.^3)/norm(x)^3
> 
>>> My idea was now to go to frequency domain (Parseval) and obtain the PSD
>>> of x^3[n] via 2-fold convolution of X[k].
>>
>> This should be valid.
> 
> See below.


Qualitatively, yes.

Here is the DFT plot of one particular instance of x[n] (drawn as above)
and x^3[n]: http://snag.gy/YFEDi.jpg

And this my "PSD" approximation: http://snag.gy/exHY9.jpg

Clearly it follows the same shape ...


>>> If I would just approximate the PSD of x[n] as ideal square function,
>>> then X[k]*X[k]*X[k] would be a quadratic and I can find an analytic
>>> estimate. Of course, this does not match at all with the numeric results
>>> when I have a random signal.
>>
>> Do you mean you've tried it and it doesn't work?  Are you sure you have 
>> enough samples?
> 
> Not sure.
> 
>> If the expression you're trying to get estimates of means what I think it 
>> does, then you're dividing by a quantity with a nonzero probability of 
>> being zero, so your expression can go to infinity -- are you sure the 
>> math isn't just blowing up in your face?
> 
> I don't think so. I don't get means but mean-squares, sum(x.^6) /
> sum(x.^2)^3. Both expressions should not be near zero (it's essentially
> "signal energy")
> 
> What I want is not to do the *numerical* estimate but derive an
> *analytical* expression for the estimate based on the assumption that
> x[n] is BLWGN described by zero mean, BW, variance.
> 
> I think my assumption is wrong to just replace my spectrum with the
> analytical PSD. On the other hand I am very sure it should be possible
> to derive an analytical expression for the estimate.
> 
> Let me quickly continue elaborating what I did:
> Given the signal generated above, I always get a similar *numerical* answer:
> 
>>> 20*log10(norm(x.^3)/norm(x)^3)
> ans =
> 
>   -74.3041
> 
> Of course, also if I do it with DFT (works due to Parseval):
> 
>>> X = fft(x);
>>> X3 = fftshift(cconv(cconv(fftshift(X),fftshift(X),N),fftshift(X),N)/N^2);
>>> 20*log10(norm(X3)/norm(X)^3*N)
> ans =
> 
>   -74.3041
> 
> Now the "idea" from above was to replace "X" with the analytical form of
> the signals PSD:
> 
> X = 0*X;
> X(1:2040) = c;
> X(18362:end) = c;
> x = ifft(X);
> 
> But now (same result in DFT domain and convolution):
> 
>>> 20*log10(norm(x.^3)/norm(x)^3)
> ans =
> 
>   -16.5779

If I use DFT domain and actually *square* X3 from above I get nearly the
same result:

X = zeros(N,1);
X(1:2040) = 1; % scaling does not matter
X(18362:end) = 1;
x = ifft(X);

X3 = (fftshift(cconv(cconv(
  fftshift(X),fftshift(X),N),fftshift(Xin),N)/N^2));

% square it
X3 = X3.^2;

>> 20*log10((norm(X3)/norm(X)^3*N))
ans =

  -73.1083

So the result is just off by 1.2 dB!

Is this just coincidence? I am aware that it is the *power* spectral
density (so squared) but I am confused why just squaring X3 but not X
for the result to match?

Also I am slightly concerned that the shapes to not match "perfectly"
(the analytical version rolls off more quickly):

http://snag.gy/rbH8T.jpg

> So I can't just do that. Probably what you mean is I would need to
> extend the signal in time somehow? Maybe there is a different approach
> at all? If it's easier to approach it in the continuous-time domain I'm
> also open for hints on that.

Regards,
Peter

On Saturday, April 4, 2015 at 12:34:57 PM UTC-7, Peter Mairhofer wrote:
...
> Mathematically I can just normalize by some period and let period go to
> infinity in the limit.
> 
> Continuing with MATLAB syntax, I would seek an (analytical) estimate for
> norm(x.^3)/norm(x)^3
> 
...
> Peter

The approach of calculating the continuous/infinite versions of signal functions with finite numbers of sets of finite numbers of samples is fundamentally flawed.

First, perhaps this just a vocabulary error, there are no periods in white Gaussian noise. There are independent blocks of samples of white Gaussian noise.

Second, the ratio of the variance of the single-block-discrete estimate of PSD to the expectation value of the discrete estimate (and the continuous/infinite value) of the PSD does not change with block size. Going to bigger n has no effect on the error.

Third, it is the finite-number-of-blocks estimate of the expectation value of the discrete PSD that converges to the expectation value of the PSD as the number of independent blocks over which the discrete PSD is calculated is increased.

Dale B. Dalrymple

On Saturday, April 4, 2015 at 1:32:02 PM UTC+5:30, Peter Mairhofer wrote:
> Hi,
> 
> I have a (for simplicity brickwall) bandlimited white Gaussian noise
> signal x[n]. Hence I know roughly its PSD.
> 
> I would like to get an estimate of ||x^3[n]||/||x[n]||^3  (||.|| l2-Norm
> of the discrete time signal x[n]).
> 
> My idea was now to go to frequency domain (Parseval) and obtain the PSD
> of x^3[n] via 2-fold convolution of X[k].
> 
> If I would just approximate the PSD of x[n] as ideal square function,
> then X[k]*X[k]*X[k] would be a quadratic and I can find an analytic
> estimate. Of course, this does not match at all with the numeric results
> when I have a random signal.
> 
> Any hints on that?
> 
> Regards,
> Peter

- x(n) is Gaussian IID (Independent - Means White, Identically distributed). E(x^2) = sigma^2.
- x(n) is filtered with a Brick Wall filter - This is "sinc" pulse in the time domain. Say the resulting sequence is y(n)
- Now y(n) is Still Gaussian and Identically distributed, but not white (Samples are no more independent).
- Since "sinc" function is absolutely integrable, Each sample y(n) is identically distributed Gaussian with variance a scale of sigma^2. The scale depends on the brick-wall filter used - Some summation of the "sinc" samples - Absolutely integrable.
- Now you want to find out sqrt(E(y^6)) which is norm(y^3) - Distribution of y is Gaussian with zero mean. E(y^p) is factorial(factorial(p-1))*(std(y))^p

On Monday, April 6, 2015 at 12:01:27 AM UTC-7, antenna wrote:
...

The original poster has postulated that he has a set of n samples: x(n) from a continuous-infinite Gaussian process x and a method of converting those n samples into n samples y(n) of a bandlimited process y.

> 
> - x(n) is Gaussian IID (Independent - Means White, Identically distributed).

No. That's the continuous-infinite x. x(n) as a finite set of samples has a discrete distribution and each time another independent block of n samples of x is generated, it will have a different discrete distribution.

> E(x^2) = sigma^2.

Now, what does this mean?
If it means :
"continuous-E(continuous-x^2) = continuous-sigma(continuous-x)^2"
that is one thing.

If it means:
"discrete-E(x(n)^2) = discrete-sigma(x(n))^2"
then please note that the value of the discrete-E is different for each independent set of n samples of x.

That means that for stochastic variables like x, "discrete-E(x(n)^2)" is an inconsistent estimator of "continuous-E(continuous-x^2)".

Further, the expected error between the continuous and discrete versions does not decrease as n becomes larger.

See:
http://mathworld.wolfram.com/ExpectationValue.html
for better notation.

> - x(n) is filtered with a Brick Wall filter - This is "sinc" pulse in the time domain. Say the resulting sequence is y(n)
> - Now y(n) is Still Gaussian and Identically distributed, but not white (Samples are no more independent).

No, y(n) is discrete and has a discrete distribution which is different each time you calculate another block of samples y(n).

> - Since "sinc" function is absolutely integrable, Each sample y(n) is identically distributed Gaussian with variance a scale of sigma^2. The scale depends on the brick-wall filter used - Some summation of the "sinc" samples - Absolutely integrable.

Absolutely integrable is of property of the theoretical continuous-infinite variable y, not of finite sets of discrete samples of y.

> - Now you want to find out sqrt(E(y^6)) which is norm(y^3) - Distribution of y is Gaussian with zero mean. E(y^p) is factorial(factorial(p-1))*(std(y))^p

Now if you want to find out "sqrt(continuous-E(continuous-y^3))" you can apply that formula to continuous-y. But if you insert finite-discrete y(n) into the right hand side of the equation, you get an inconsistent estimator of "sqrt(continuous-E(continuous-y^3))".

If you want to calculate a consistent estimator of a continuous property of a continuous infinite stochastic variable, you cannot do it with a single finite set of samples from the stochastic variable. But you might be able to approximate it from a large number of blocks of finite numbers of discrete samples.

Discussions of calculations with samples from stochastic variables can get drearily verbose, but if you drop the labels, someone always seems to make a continuous-infinite assumption on one side of their equation and put finite -discrete values on the other, then ask "What's wrong with this?".

Dale B Dalrymple

On Sat, 04 Apr 2015 12:34:54 -0700, Peter Mairhofer wrote:

> On 2015-04-04 7:57, Tim Wescott wrote:
>> On Sat, 04 Apr 2015 01:01:59 -0700, Peter Mairhofer wrote:
>> 
>>> Hi,
>>>
>>> I have a (for simplicity brickwall) bandlimited white Gaussian noise
>>> signal x[n]. Hence I know roughly its PSD.
>>>
>>> I would like to get an estimate of ||x^3[n]||/||x[n]||^3  (||.||
>>> l2-Norm of the discrete time signal x[n]).
>> 
>> Your notation isn't clear here.  Are you referring to x[n] as the n'th
>> sample of x, here, or are you referring to x[n] as the whole signal? 
>> If it's the "whole signal" and the signal is of infinite extent, I
>> don't think the answer is finite.
> 
> It is the whole signal.
> 
> For simplicity I assume now I have a signal of, say, length N=20400 with
> the above properties.
> 
> One instance of the signal is generated with:
> 
> % WGN x = sqrt(sigma2) * randn(N,1);
> 
> % BLWGN, brickwall filter X = fft(x);
> X(2041:18361) = 0;
> x = ifft(X);
> 
> Mathematically I can just normalize by some period and let period go to
> infinity in the limit.
> 
> Continuing with MATLAB syntax, I would seek an (analytical) estimate for
> norm(x.^3)/norm(x)^3
> 
>>> My idea was now to go to frequency domain (Parseval) and obtain the
>>> PSD of x^3[n] via 2-fold convolution of X[k].
>> 
>> This should be valid.
> 
> See below.
> 
>>> If I would just approximate the PSD of x[n] as ideal square function,
>>> then X[k]*X[k]*X[k] would be a quadratic and I can find an analytic
>>> estimate. Of course, this does not match at all with the numeric
>>> results when I have a random signal.
>> 
>> Do you mean you've tried it and it doesn't work?  Are you sure you have
>> enough samples?
> 
> Not sure.
> 
>> If the expression you're trying to get estimates of means what I think
>> it does, then you're dividing by a quantity with a nonzero probability
>> of being zero, so your expression can go to infinity -- are you sure
>> the math isn't just blowing up in your face?
> 
> I don't think so. I don't get means but mean-squares, sum(x.^6) /
> sum(x.^2)^3. Both expressions should not be near zero (it's essentially
> "signal energy")

Blegh.  I read your equation, and then misinterpreted it.  That's entirely 
my fault.

Yes, you should get something sensible.  Or at least if you normalize 
correctly for the vector length you should.

> What I want is not to do the *numerical* estimate but derive an
> *analytical* expression for the estimate based on the assumption that
> x[n] is BLWGN described by zero mean, BW, variance.
> 
> I think my assumption is wrong to just replace my spectrum with the
> analytical PSD. On the other hand I am very sure it should be possible
> to derive an analytical expression for the estimate.

It should be quite possible to do this analytically.  It may even be 
tractable if you have the math chops for it.  I'm not sure what your 
background is, though -- you're talking graduate-level statistics and 
stochastic systems here.

If you know that x[n] is a filtered version of a zero-mean Gaussian random 
signal, and if you know the filter that generates x[n], then you know 
everything you need:

- You know the pdf of x_n for any given n
- you know the correlation of x_n and x_{n-m} for any given n and m

From that, you can calculate the pdf of the numerator and denominator of 
your metric, which means that you can, in turn, calculate the pdf of your 
metric, and then you have its mean, variance, kurtosis, possibly even its 
halitosis, and every other statistic that you may wish to have on it.

This will be easiest to do with discrete-time signals and a FIR filter, 
because the correlation of x_n with x_{n-m} will be zero for m larger than 
your FIR length.  That means you'll be dealing with finite sums.  Next 
harder will be IIR filters and discrete-time signals (infinite sums), and 
after that will be continuous-time signals, where you'll get to experience 
the delights of solving integral equations.

I would find the continuous-time signal case a delight, but I'd only do it 
if I had lots and lots of time on my hands, if someone said "yes, I know 
it's probably not immediately necessary, but I'll pay for it", or if I 
were a university prof (in which case the statement about immediate 
necessity and payment would be implied).

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

On 2015-04-06 0:01, antenna wrote:
> On Saturday, April 4, 2015 at 1:32:02 PM UTC+5:30, Peter Mairhofer wrote:
>> Hi,
>>
>> I have a (for simplicity brickwall) bandlimited white Gaussian noise
>> signal x[n]. Hence I know roughly its PSD.
>>
>> I would like to get an estimate of ||x^3[n]||/||x[n]||^3  (||.|| l2-Norm
>> of the discrete time signal x[n]).
>>
>> My idea was now to go to frequency domain (Parseval) and obtain the PSD
>> of x^3[n] via 2-fold convolution of X[k].
>>
>> If I would just approximate the PSD of x[n] as ideal square function,
>> then X[k]*X[k]*X[k] would be a quadratic and I can find an analytic
>> estimate. Of course, this does not match at all with the numeric results
>> when I have a random signal.
>>
>> Any hints on that?
>>
>> Regards,
>> Peter
> 
> - x(n) is Gaussian IID (Independent - Means White, Identically distributed). E(x^2) = sigma^2.
> - x(n) is filtered with a Brick Wall filter - This is "sinc" pulse in the time domain. Say the resulting sequence is y(n)
> - Now y(n) is Still Gaussian and Identically distributed, but not white (Samples are no more independent).
> - Since "sinc" function is absolutely integrable, Each sample y(n) is identically distributed Gaussian with variance a scale of sigma^2. The scale depends on the brick-wall filter used - Some summation of the "sinc" samples - Absolutely integrable.
> - Now you want to find out sqrt(E(y^6)) which is norm(y^3) - Distribution of y is Gaussian with zero mean. E(y^p) is factorial(factorial(p-1))*(std(y))^p

Thank you! That was great and helpful!

Going through time domain is actually much simpler than frequency
domain/PSD.

Interesting that for ||x^3||/(sigma^2 ||x||) it is bandwidth independent
(given that variance of y is known) and the variance cancels out in the
end (anyway) - so it is always a constant: 15 (matching with my
empirical results).

Peter