| Re: Wiener filter [message #28528 is a reply to message #28424] |
Wed, 19 December 2001 13:48   |
James Kuyper Jr.
Messages: 10
Registered: November 2001
|
Junior Member |
|
|
Richard Tyc wrote:
>> Optimal Wiener filtering of a one-dimensional data set is described in
>> section 12.6 of "Numerical Recipes in C", by Preuss et.al. It cites
>> three books on signal processing as references. The basic result is that
>> if you have a corrupted signal with the fourier spectrum S(f),
>> containing noise with a fourier specturm N(f), it can be shown
>> rigorously that the optimal (in the sense of a least-squares fit)
>> frequency filter for removing the noise is:
>>
>> |S(f)|^2
>> phi(f) = --------------------
>> |S(f)|^2 + |N(f)|^2
>>
>
>
> What exactly is |S(f)|^2
It's proportional to the power spectral density of the corrupted signal.
S(f) is the fourier transform of the the corrupted signal.
> If I have a 2D corrupted image, say I(x,y)
>
> Is it ABS( FFT(I) ) ^ 2 or the magnitude of the complex FFT result
> squared (Power Spectrum) squared ?
It's ABS(FFT(I))^2. I think that the following code should be a more
efficient way to calculate the same value:
s = FFT(I)
ps = DOUBLE(s,conj(s))
...
> It seems some knowledge of the noise is required. ...
Correct.
> ... What if it was modeled as
> 'white noise' where it would be constant at all spatial frequencies.
Optimal Wiener filtering is a way of using your knowledge of the
frequency power spectrum of the noise, to extract it from the data. If
you're using a "white noise" spectrum because you know that your noise
has that characteristic, that's reasonable. If you're using "white
noise" because you're not sure what the noise spectrum looks like, and
are afraid to commit yourself, Wiener filtering is inappropriate. Keep
in mind that you need to know not merely the frequency dependence of the
noise, but also it's absolute magnitude.
On the other hand, you don't need to know the noise power spectrum very
precisely. The result of the filtering is insensitive to small errors in
the assumed noise spectrum, just like sqrt(x) is insensitve to small
errors in 'x'.
> A paper I am using that discusses this in the context of my problem points
> out, "....Assuming that noise power spectrum is white, the mean spectral
> density at high spatial frequencies was calculated and subtracted from P(f)
> (the power spectral density of the corrupted image) to estimate S(f) (power
> spectral density of uncorrupted image). Can you shed any light on this in
> terms of IDL code ??
Note a difference in notation here: that quotation identifies S(f) as
the power spectral density of the uncorrupted image. I was using that
same notation to mean the fourier transform of the corrupted signal.
I'm afraid I can't convert that to code; the only tricky step is one
that they've given no details about. That step is the one where they
estimate the power spectrum of the noise. They said that they were
assuming a "white noise" spectrum, but that still leaves them with the
problem of estimating the amplitude of the noise. One plausible approach
is to plot the power spectrum of your signal, and decide to model it as
the sum of two simple curves with known shapes. Then use regress() to
fit the data to a linear sum of those two curves.
The part they do explain is trivial. Using the notation from that quote,
P(f) = S(f) - N(f), where P(f) is the power spectrum of the uncorrupted
signal, S(f) is the power specttrum of the corrupted signal, and N(f) is
the power spectrum of the noise. (Note change of notation from previous
context).
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
