Inquiry about white noise and wiener deconvolutions

Started by February 24, 2009
I&rsquo;m interested in trying a wiener deconvolution for some
electro-physiology data that I have recorded, and I&rsquo;m not sure how to
estimate the white noise power for use in that. According to Wikipedia, the
proper terms to use are the PSD of both the noise and the true signal in
the equation used for the wiener deconvolution. (As an aside question, I
was going to use the recorded signal as a proxy for the true signal. Is
that common or are there better ways of going about this?) The analog
filters completely cutoff the highest frequency portion of the spectrum
prior to digitization, in a manner such that there is a range of high
frequencies below the Nyquist in which the signal is completely absent. I
was planning on using the spectral power averaged over those frequencies as
an estimate of the noise.

However- the way I envision the recording system working is that following
the application of the analog filters, something akin to a zero-mean
Gaussian random value with some particular variance is added to each
recorded sample, with the amount added to any given sample being
independent of the amount added to any other given sample. In matlab, this
would be akin to something like:
Y= X + sigma*randn( length(X) ,1 )
Where X is the filtered signal, and Y is the recorded signal. Although the
PSD of the noise, i.e. randn( length(X) ,1 ), is flat with respect to
frequency, the value at any particular frequency depends on the length of
X, growing smaller as the length of X grows. A simple matlab simulation
shows this. My take on this is that the randn variable has a certain finite
constant power based on its overall variance, so that when the signal goes
on for more and more time, that finite variances get spread more and more
thinly across more discrete frequencies.

So then it seems that when the wiener filter is applied, the PSD of such
noise relative to the signal would grow smaller and smaller, rendering the
wiener filter eventually useless for a long enough signal. Note that
assuming the signal is always present and not horrifically non-stationary,
it&rsquo;s PSD will be independent of it&rsquo;s length, unlike the PSD of the
noise. Should the mean-square spectrum of the noise and signal be used
instead? Or am I getting something wrong? Please let me know if anything I
said above is incorrect.

On Feb 24, 7:27&#2013266080;pm, "john917" <nelsonmj...@hotmail.com> wrote:
> I&#2013266066;m interested in trying a wiener deconvolution for some > electro-physiology data that I have recorded, and I&#2013266066;m not sure how to > estimate the white noise power for use in that. According to Wikipedia, the > proper terms to use are the PSD of both the noise and the true signal in > the equation used for the wiener deconvolution. (As an aside question, I > was going to use the recorded signal as a proxy for the true signal. Is > that common or are there better ways of going about this?) The analog > filters completely cutoff the highest frequency portion of the spectrum > prior to digitization, in a manner such that there is a range of high > frequencies below the Nyquist in which the signal is completely absent. I > was planning on using the spectral power averaged over those frequencies as > an estimate of the noise. > > However- the way I envision the recording system working is that following > the application of the analog filters, something akin to a zero-mean > Gaussian random value with some particular variance is added to each > recorded sample, with the amount added to any given sample being > independent of the amount added to any other given sample. In matlab, this > would be akin to something like: > &#2013266080; &#2013266080; &#2013266080; &#2013266080; Y= X + sigma*randn( length(X) ,1 ) > Where X is the filtered signal, and Y is the recorded signal. Although the > PSD of the noise, i.e. randn( length(X) ,1 ), is flat with respect to > frequency, the value at any particular frequency depends on the length of > X, growing smaller as the length of X grows. A simple matlab simulation > shows this. My take on this is that the randn variable has a certain finite > constant power based on its overall variance, so that when the signal goes > on for more and more time, that finite variances get spread more and more > thinly across more discrete frequencies. > > So then it seems that when the wiener filter is applied, the PSD of such > noise relative to the signal would grow smaller and smaller, rendering the > wiener filter eventually useless for a long enough signal. Note that > assuming the signal is always present and not horrifically non-stationary, > it&#2013266066;s PSD will be independent of it&#2013266066;s length, unlike the PSD of the > noise. Should the mean-square spectrum of the noise and signal be used > instead? Or am I getting something wrong? Please let me know if anything I > said above is incorrect.
Your observation with respect to the noise power per bin decreasing as you increase the length of your signal (and therefore the length of FFT that you're using to estimate the PSD) is correct. Think about it this way: your AWGN has a finite average power (variance). That power is uniformly distributed across all frequencies. The FFT is just a filterbank that conveniently places evenly-spaced filters across each band. If the noise has power P, then by using an FFT of length 16 to estimate the PSD, you would expect to see power levels of P/16 in each bin. If you use a longer FFT though, say 64, then you would expect a power level of P/64 in each bin. The noise still has the same total power; you're just splitting it into smaller sections when you use a longer FFT for the PSD estimate. In this way, it's helpful to think of the filterbank interpretation of the FFT that I suggested above; the value of each FFT bin really is just a collection of critically-decimated outputs of discrete-time boxcar filters whose inputs have been downconverted from their respective bin center frequencies. You can then use those outputs to estimate the power level in each of the subband outputs. Jason
On Tue, 24 Feb 2009 18:27:16 -0600, john917 wrote:

> I&rsquo;m interested in trying a wiener deconvolution for some > electro-physiology data that I have recorded, and I&rsquo;m not sure how to > estimate the white noise power for use in that. According to Wikipedia, > the proper terms to use are the PSD of both the noise and the true > signal in the equation used for the wiener deconvolution. (As an aside > question, I was going to use the recorded signal as a proxy for the true > signal. Is that common or are there better ways of going about this?)
If you don't know much about the processes involved that's not a bad way to go. You may know all the numbers for other reasons -- the signal PSD is well known, the noise power comes from a known process, etc. But deducing your parameters from measurements is perfectly valid.
> The analog filters completely cutoff the highest frequency portion of > the spectrum prior to digitization, in a manner such that there is a > range of high frequencies below the Nyquist in which the signal is > completely absent. I was planning on using the spectral power averaged > over those frequencies as an estimate of the noise.
This sounds reasonable.
> > However- the way I envision the recording system working is that > following the application of the analog filters, something akin to a > zero-mean Gaussian random value with some particular variance is added > to each recorded sample, with the amount added to any given sample being > independent of the amount added to any other given sample.
You are describing zero mean Gaussian white noise.
> In matlab, > this would be akin to something like: > Y= X + sigma*randn( length(X) ,1 ) > Where X is the filtered signal, and Y is the recorded signal. Although > the PSD of the noise, i.e. randn( length(X) ,1 ), is flat with respect > to frequency, the value at any particular frequency depends on the > length of X, growing smaller as the length of X grows. A simple matlab > simulation shows this. My take on this is that the randn variable has a > certain finite constant power based on its overall variance, so that > when the signal goes on for more and more time, that finite variances > get spread more and more thinly across more discrete frequencies.
I suspect that you are conflating Matlab with reality, and confusing a PSD with the result of an FFT on a sample of a noise process. The notion of power spectral _density_ is that it puts so many watts (or so many signal-units-squared) into each unit of bandwidth -- hence, it is expressed as watts/Hz, or by its square root as (whatever natural units)/ rt-Hz. How your FFT resolves this noisy signal depends entirely on how your math package scales its FFT output. Fortunately for you, every writer of FFT algorithms on earth scales their outputs in the very best possible way. Unfortunately for you, each one of them feels that "the best" is something different from the next guy. So what Matlab does with your signal is nearly immaterial, until you go to interpret your PSD results. Your best bet is to find the total energy of your sample vector in the time domain, find the sum-of-squares of the energy of the FFT of the sample vector, then use Parseval's theorem to figure out the scaling factor that squares everything up.
> > So then it seems that when the wiener filter is applied, the PSD of such > noise relative to the signal would grow smaller and smaller, rendering > the wiener filter eventually useless for a long enough signal. Note that > assuming the signal is always present and not horrifically > non-stationary, it&rsquo;s PSD will be independent of it&rsquo;s length, unlike the > PSD of the noise. Should the mean-square spectrum of the noise and > signal be used instead? Or am I getting something wrong? Please let me > know if anything I said above is incorrect.
I think you are getting something wrong -- see my discussion above. Logically, the way to verify this is to consider a recording of no signal at all. In that case your Matlab model of the recording will be Y= X + sigma*randn( length(X) ,1 ) with X = 0, or just Y= sigma*randn( length(X) ,1 ) You can tell just by looking at this that the expected sum of squares of the samples of Y (i.e. the energy in Y) will increase linearly with the length of the sample. -- http://www.wescottdesign.com