Whitening in Frequency Domain

Hello - 

I want to create a whitening filter in the frequency domain, and I have
implemented such a filter (as an adaptation of some older code), but I
believe it is incorrect as I am getting some odd edge effects.  I have come
here to hope to receive some help on my theoretical understanding of the
process.

Basically, I want to create the following filter H = 1/sqrt(|X|) in the
frequency domain, so I perform the following:

1. take first N samples [1:N] of the input signal (x[n])
2. take fft to obtain X
3. calculate H = 1/sqrt(|X|)
4. multiply in freq domain X_w = X*H
5. take inverse fft of X_w to get whitened time domain x_w[n].
6. keep first N/2 samples of x_w[n] as output
7. take x[N/2:N/2+N] and repeat 2-6

My issue has to do with the multiplication in the frequency domain. 
Multiplication in the frequency domain is performing circular convolution,
and the overlap-save technique tells us that if the input frame is size N
and the filter length is size L, then L-1 samples are invalid due to time
aliasing, leaving N - (L-1) valid samples.

Now, in the filter that I have built, I believe the filter length AND the
input frame are both N, thus resulting in N-N+1=1 valid output sample?

My first thought was to make H of length N/2, but I'm not quite sure how to
do that.  Basically, my thought was that I need to ensure that h has N/2
nonzero values and N/2 zero values at the end.  However, I am not sure how
to accomplish this.

Clear as mud?

I guess my question is - How do I create a whitening filter in frequency
domain that is free of time aliasing.

Thanks.

On Thursday, March 7, 2013 7:39:30 PM UTC-8, dspillini wrote:
> Hello - 
> ...
> I guess my question is - How do I create a whitening filter in frequency
> domain that is free of time aliasing.
> 
> Thanks.

First learn about frequency domain filtering and pay attention to "zero fill". Consider the whole chapter starting at:
https://ccrma.stanford.edu/~jos/sasp/Overlap_Add_OLA_STFT_Processing.html
and
http://www-dsp.elet.polimi.it/ispg/images/pdf/audio/materiale/ola.pdf

Then figure out what form of smoothed spectral estimate your application field prefers for whitening.

Dale B. Dalrymple

>I guess my question is - How do I create a whitening filter in frequency
>domain that is free of time aliasing.

FFT-iFFT is not magic, it is just subband filters with
decimation-interpolation:
https://ccrma.stanford.edu/~jos/sasp/Downsampled_STFT_Filter_Bank.html
https://ccrma.stanford.edu/~jos/sasp/Filter_Bank_Reconstruction.html
To avoid aliasing just use proper FIR prototype(window(should be several
times larger than the block FFT -> polyphase filter)) and
decimation-interpolation coefficient(overlapping blocks).

Look for example simulink model
http://electronix.ru/forum/index.php?=&showtopic=23652&view=findpost&p=929325

You can't do bin-by-bin frequency response manipulation like you are proposing without having aliasing effects. You will need to use a smoothed version of H.

One way to do this is to calculate H as you are presently doing, then inverse transform back to the time domain, and the apply a window function where the length of the window is N/2.  Then forward transform back to the frequency domain and use this new H function to multiply your frames. The new H function will be a smoothed version of the original.

Full disclosure, I have never tried this myself so I'm not sure if my theory is correct, but it has a certain feel of truthiness to it :)

Bob

All - thanks for the suggestions.  I'll do some more research re: your
resources and the STFT and get reply back with anymore questions I
encounter.

On Friday, March 8, 2013 4:08:22 AM UTC-8, radam...@gmail.com wrote:
> You can't do bin-by-bin frequency response manipulation like you are proposing without having aliasing effects. You will need to use a smoothed version of H.

Yes

> One way to do this is to calculate H as you are presently doing, then inverse transform back to the time domain, and the apply a window function where the length of the window is N/2.  Then forward transform back to the frequency domain and use this new H function to multiply your frames. The new H function will be a smoothed version of the original.
> 
> Full disclosure, I have never tried this myself so I'm not sure if my theory is correct, but it has a certain feel of truthiness to it :)
> 
> Bob

Well, any of the cosine-sum windows (Hamming, von Hann, Blackman, Blackman-harris, any of the windows in the Nuttall paper) can be applied as a small convolution in the frequency domain. For a N term window it's a 2N-1 term convolution, so 3 point convolution for Hamming or von Hann. This is much less processing than the ifft/fft processing and windowing but exactly equivalent. But once people have to decide what coefficients they need to use, they often just use a frequency domain sliding boxcar to minimize multiplies or to remove bias in the estimate of strong components: a split boxcar, or to remove arbitrarily strong components from a background estimate: a sliding median.

Dale B. Dalrymple

On Thu, 07 Mar 2013 21:39:30 -0600, dspillini wrote:

> Hello -
> 
> I want to create a whitening filter in the frequency domain, and I have
> implemented such a filter (as an adaptation of some older code), but I
> believe it is incorrect as I am getting some odd edge effects.  I have
> come here to hope to receive some help on my theoretical understanding
> of the process.
> 
> Basically, I want to create the following filter H = 1/sqrt(|X|) in the
> frequency domain, so I perform the following:
> 
> 1. take first N samples [1:N] of the input signal (x[n]) 2. take fft to
> obtain X
> 3. calculate H = 1/sqrt(|X|)
> 4. multiply in freq domain X_w = X*H
> 5. take inverse fft of X_w to get whitened time domain x_w[n]. 6. keep
> first N/2 samples of x_w[n] as output 7. take x[N/2:N/2+N] and repeat
> 2-6
> 
> My issue has to do with the multiplication in the frequency domain.
> Multiplication in the frequency domain is performing circular
> convolution, and the overlap-save technique tells us that if the input
> frame is size N and the filter length is size L, then L-1 samples are
> invalid due to time aliasing, leaving N - (L-1) valid samples.
> 
> Now, in the filter that I have built, I believe the filter length AND
> the input frame are both N, thus resulting in N-N+1=1 valid output
> sample?
> 
> My first thought was to make H of length N/2, but I'm not quite sure how
> to do that.  Basically, my thought was that I need to ensure that h has
> N/2 nonzero values and N/2 zero values at the end.  However, I am not
> sure how to accomplish this.
> 
> Clear as mud?
> 
> I guess my question is - How do I create a whitening filter in frequency
> domain that is free of time aliasing.

This sounds like yet another attempt to use the FFT as a magic wand, 
hoping to get the incantations right.

What are you _really_ trying to do?  Given that you're talking about 
whitening, I assume that somewhere buried underneath whatever razz-matazz 
you're dealing with is a desire to make a Wiener filter.

Is that true?

Why, if you're doing signal estimation, do you feel that you must use an 
FFT?

If you don't feel that you must use an FFT, why do you feel that an FFT 
is best?

-- 
My liberal friends think I'm a conservative kook.
My conservative friends think I'm a liberal kook.
Why am I not happy that they have found common ground?

Tim Wescott, Communications, Control, Circuits & Software
http://www.wescottdesign.com

>You can't do bin-by-bin frequency response manipulation like you are
propos=
>ing without having aliasing effects. You will need to use a smoothed
versio=
>n of H.
>
>One way to do this is to calculate H as you are presently doing, then
inver=
>se transform back to the time domain, and the apply a window function
where=
> the length of the window is N/2.  Then forward transform back to the
frequ=
>ency domain and use this new H function to multiply your frames. The new H
=
>function will be a smoothed version of the original.
>
>Full disclosure, I have never tried this myself so I'm not sure if my
theor=
>y is correct, but it has a certain feel of truthiness to it :)
>
>Bob
>
>

Bob - 

I had this thought as well, however how do I guarantee that this way would
provide me with constant group delay?  My original H is purely real with
phase = 0, thus a group delay of 0.  With what your proposing, a group
delay of zero, or even a constant group delay would not be guaranteed,
correct?

Thanks.

On 3/8/13 1:56 PM, dspillini wrote:
>> You can't do bin-by-bin frequency response manipulation like you are propos=
>> ing without having aliasing effects. You will need to use a smoothed versio=
>> n of H.
>>
>> One way to do this is to calculate H as you are presently doing, then inver=
>> se transform back to the time domain, and the apply a window function where=
>> the length of the window is N/2.  Then forward transform back to the frequ=
>> ency domain and use this new H function to multiply your frames. The new H
> =
>> function will be a smoothed version of the original.
>>
>> Full disclosure, I have never tried this myself so I'm not sure if my theor=
>> y is correct, but it has a certain feel of truthiness to it :)
>>

Bob, i wish you might find out what it is in your email client that does 
not put in nice line breaks at all or not *between* words.  i end up 
having to fiddle with your quoted text all the time.


>
> I had this thought as well, however how do I guarantee that this way would
> provide me with constant group delay?  My original H is purely real with
> phase = 0, thus a group delay of 0.  With what your proposing, a group
> delay of zero, or even a constant group delay would not be guaranteed,
> correct?

you can make a constant group delay (which is the same as a constant 
phase delay where the constant delay is the same for both).

whatever you have in the frequency domain will be Hermitian symmetric, 
right?  and since your original H is "zero phase", then when you inverse 
DFT, that result will be both real and symmetric about n=0.  put the 
latter half ( N/2 < n < N) in front of the first half (because it's like 
"negative time" anyway and you should have even symmetry. apply a 
even-symmetric window of your choice to shorten it.  then to make the 
implied FIR of your fast-convolver causal, delay it by enough so that 
there is nothing non-zero in the "negative time" section, then swap it 
back (this is what fftshift() in MATLAB does) and DFT it for your 
"smoothed" frequency response.

you will have a causal FIR filter, that is linear phase (constant group 
and phase delay), causal, with enough zeros padded in it so that you can 
accomplish linear convolution by use of circular convolution.

it need not be N/2 in length, but the shorter you can make the FIR the 
longer you can make the hop length and the more samples you can process 
per frame.

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

>What are you _really_ trying to do?  Given that you're talking about 
>whitening, I assume that somewhere buried underneath whatever razz-matazz

>you're dealing with is a desire to make a Wiener filter.

I am trying to remove interference from my signal of interest.  The SOI is
bursty, so there are periods where it is off.  So, during those OFF
periods, I want to estimate the spectra in an attempt to "whiten" the
interferers.  In order to do so, I am performing the following:
1. set ON periods to 0 to obtain x_gate
2. take fft of x_gate to get X_gate
3. take fft of x to get X
4. calculate H=1/|X_gate|
5. multiply X*H
6. perform IFFT to get y = IFFT(X*H)

>Why, if you're doing signal estimation, do you feel that you must use an 
>FFT?
>
>If you don't feel that you must use an FFT, why do you feel that an FFT 
>is best?

I was only using an FFT so I could estimate the spectra in order to obtain
the inverse magnitude (1/|X_gate|).