Forums

prewhitening a signal in matlab

Started by mdma...@gmail.com October 6, 2010
Hello,
I am looking at a time series that doesn't have a white spectrum, and I would like to whiten it. In theory, I understand whitening as taking a signal which is not white, fourier transforming it, then dividing its magnitude by its magnitude (so that it is 1 = white) and then inverse fourier transforming it.
I know that for a signal which is discrete and finite, that this division by its magnitude spectrum doesn't make sense - although I cannot justify the real reasons why it doesn't make sense (perhaps you may help me in understanding this too). Instead then, one should compute a "smoothed" magnitude and divide the original magnitude by the "smoothed" magnitude and get a final magnitude that is approximately equal to 1 but not constant for all frequencies.
So far so good. I was using matlab's pmtm function to smooth in a program that looked like the following:
ptot=dlmread(‘time_series.txt’);
dt=0.001;
N=length(ptot);
df=1 / (1*(dt)*(N-1)); % Hz
p_fft_abss(fft(ptot,N));
window=hann(N);
ptot=ptot.*window;

%%%% SMOOTHING SECTION WITH PMTM (multitaper Thomson algorithm)
[Pxx,w] = pmtm(ptot,3,length(ptot));
Pxx=Pxx';
w=w';
wHz=w./(2*pi*dt); % this takes us in Hz
% construct other half-sides
Pxx_left=fliplr(Pxx(2:length(Pxx)));
Pxx_total=[Pyy_left Pxx];
wHz_left=fliplr(wHz(2:length(wHz)));
wHz_total=[-wHz_left wHz];

Pxx_phase=angle(fft(ptot,N)); % note: phase doesn't change with smoothing
Pxx_mag=sqrt(ifftshift(Pxx_total)); % THIS IS THE SMOOTHED magnitude

%%% SCALE THE ORIGINAL AND SMOOTHED AMPLITUDES
p_fft_abs=p_fft_abs./max(p_fft_abs);
Pxx_mag=Pxx_mag./max(Pxx_mag);

%%% WHITENING SECTION: we divide the original magnitude with the smoothed one
white_mag=p_fft_abs./(Pxx_mag);
white_dom_freq=white_mag.*exp(i*Pxx_phase); % this is the 'prewhitened' signal in frequency space
white_dom_time=ifft(white_dom_freq); %this is the prewhitened signal in time domain

However, I noticed that the resulting “white” signal in time that I had created had certain “jumps” at its beginning and end in the time domain and I don’t really know why. I am hoping you can cross-check if my code should do what I want it to do. A further question I have concerns the windowing part of the signal: I understand why windowing in time domain is important to limit spectral leakage in the frequency domain, but should I also be windowing the frequency domain prior to ifft in order to avoid a convolution by sinc in the time domain? Most entries I found on windowing and the ifft refer to compensating the effect of the windowing in the time domain (the window applied to the time series prior to fft), but do not refer to what a rectangular window in the frequency does to the time domain (so a window applied to the frequency domain prior to ifft).

Thank you for your support!!
Maria-Daphne
Maria-Daphne-

> I am looking at a time series that doesn't have a white spectrum, and I would like to whiten it. In theory, I
> understand whitening as taking a signal which is not white, fourier transforming it, then dividing its magnitude by
> its magnitude (so that it is 1 = white) and then inverse fourier transforming it.
> I know that for a signal which is discrete and finite, that this division by its magnitude spectrum doesn't make sense
> - although I cannot justify the real reasons why it doesn't make sense (perhaps you may help me in understanding this
> too). Instead then, one should compute a "smoothed" magnitude and divide the original magnitude by the "smoothed"
> magnitude and get a final magnitude that is approximately equal to 1 but not constant for all frequencies.
> So far so good. I was using matlab's pmtm function to smooth in a program that looked like the following:

If you make the magnitude perfectly flat, then after inverse FFT you would have an infinite impulse (Dirac Delta
function) in the time domain.

I suggest that the objective when whitening a signal would be to add energy at frequencies *other* than those carrying
the signal's basic content, and such process would not result in a perfectly flat magnitude response.

-Jeff

> ptot=dlmread(‘time_series.txt’);
> dt=0.001;
> N=length(ptot);
> df=1 / (1*(dt)*(N-1)); % Hz
> p_fft_abss(fft(ptot,N));
> window=hann(N);
> ptot=ptot.*window;
>
> %%%% SMOOTHING SECTION WITH PMTM (multitaper Thomson algorithm)
> [Pxx,w] = pmtm(ptot,3,length(ptot));
> Pxx=Pxx';
> w=w';
> wHz=w./(2*pi*dt); % this takes us in Hz
> % construct other half-sides
> Pxx_left=fliplr(Pxx(2:length(Pxx)));
> Pxx_total=[Pyy_left Pxx];
> wHz_left=fliplr(wHz(2:length(wHz)));
> wHz_total=[-wHz_left wHz];
>
> Pxx_phase=angle(fft(ptot,N)); % note: phase doesn't change with smoothing
> Pxx_mag=sqrt(ifftshift(Pxx_total)); % THIS IS THE SMOOTHED magnitude
>
> %%% SCALE THE ORIGINAL AND SMOOTHED AMPLITUDES
> p_fft_abs=p_fft_abs./max(p_fft_abs);
> Pxx_mag=Pxx_mag./max(Pxx_mag);
>
> %%% WHITENING SECTION: we divide the original magnitude with the smoothed one
> white_mag=p_fft_abs./(Pxx_mag);
> white_dom_freq=white_mag.*exp(i*Pxx_phase); % this is the 'prewhitened' signal in frequency space
> white_dom_time=ifft(white_dom_freq); %this is the prewhitened signal in time domain
>
> However, I noticed that the resulting “white” signal in time that I had created had certain “jumps” at its
> beginning and end in the time domain and I don’t really know why. I am hoping you can cross-check if my code should
> do what I want it to do. A further question I have concerns the windowing part of the signal: I understand why
> windowing in time domain is important to limit spectral leakage in the frequency domain, but should I also be
> windowing the frequency domain prior to ifft in order to avoid a convolution by sinc in the time domain? Most entries
> I found on windowing and the ifft refer to compensating the effect of the windowing in the time domain (the window
> applied to the time series prior to fft), but do not refer to what a rectangular window in the frequency does to the
> time domain (so a window applied to the frequency domain prior to ifft).
>
> Thank you for your support!!
> Maria-Daphne
Maria-Daphne-

> Hello Jeff,
>
> Thank you very much for your help. Below are some of my comments:
>
>>If you make the magnitude perfectly flat, then after inverse FFT you would
> have an >infinite impulse (Dirac Delta
>>function) in the time domain.
>
> That wouldnt be the case however, if you had a random phase right? So if
> my original signal is random, and I only make the magnitude spectrum 1, but
> dont change the phase, I shouldnt be getting an impulse correct?

Yes you could also have a flat spectrum due to random values (white noise) in the time domain. But again, if you
actually did that your original signal content would be destroyed.

-Jeff
Hello Michael,
Thank you very much for your help! Below are my comments:

>What is your application? Are you trying to model this signal using ARIMA-type >models?
My application is looking at geophone recordings of ambient noise vibrations. In seismology, some theory for recovering the earth’s structure relies on cross-correlating these ambient noise recordings. This theory relies on the assumption that these noise recordings are for earthquakes that are random in space and time, ie that the time series should be stationary. This is not the case in what we record however. Moreover, the recorded time series often has peaks at specific frequencies which overwhelm the rest of the frequencies. What is then common to seismologists is to take this signal prior to doing the cross-correlations and applying some type of spectral normalization. This will in turn broaden the band of noise signal in the cross-correlation and combat degradation caused by monochromatic persistent sources.

>If you normalize in the manner you suggest, on a frequency-by-frequncy basis, you >indeed remove the relative frequency relationships. This does whiten the sequence, >but any intermediate information is lost. Could you not simply simulate a white->noise signal directly?
This spectral normalization is commonly accomplished– within the seismological community I mean – by prewhitening the signal. I have in this context seen two prewhitening schemes: one is to add a constant over the entire spectrum of the signal (similar to what Jeff suggested (“I suggest that the objective when whitening a signal would be to add energy at frequencies *other* than those carrying
the signal's basic content, and such process would not result in a perfectly flat magnitude response.”) and the other is to take the magnitude spectrum, and divide it by itself (to get one). However, most researchers never do this substitution of the spectrum by a constant, but instead divide it by a smoothed magnitude (what I have tried to do in my code). I suspect they don’t substitute by a constant because that would be changing the data too much, ie one would loose too much information by changing completely the spectrum instead of ‘quieting’ the undesired overshoots in a frequency spectrum. Or perhaps they don’t do it because of another reason (again I suspect?): given that the signal is bandlimited if you substitute its spectrum with a constant that is like having a boxcar in the frequency domain, so the signal perhaps would be smeared in time??
Which brings me to another question: as a student I remember being advised to window my data (by blackmann hanning etc) prior to fft so to increase spectral resolution (which I see as: if you don’t window in the time domain, you are multiplying by a boxcar, which is convolving with a sinc in the frequency domain). However, why wasn’t advised to window (by blackmann hanning etc) in the frequency domain prior to ifft to improve time resolution? Encorporated in this problem of prewhitening, shouldn’t I be dividing by the smoothed spectrum, then multiplying the magnitude with a Blackman and then ifft? And leave the phase intact the entire time right?
As far as ARMA modeling I have seen in seismology applications in which these models are used to simulate earthquake ground motions, however not for the purpose of making a geophone record “white”. I wonder why they are not used in that context… I will read further into the literature you proposed. Thank you again for your help! Maria-Daphne

Hello Jeff,
Thank you very much for your help. Below are some of my comments:
>If you make the magnitude perfectly flat, then after inverse FFT you would have an >infinite impulse (Dirac Delta
>function) in the time domain.
That wouldn’t be the case however, if you had a random phase right? So if my original signal is random, and I only make the magnitude spectrum 1, but don’t change the phase, I shouldn’t be getting an impulse correct?

>I suggest that the objective when whitening a signal would be to add energy at >frequencies *other* than those carrying
>the signal's basic content, and such process would not result in a perfectly flat >magnitude response.
Yes, I think this is right, that would be equivalent to adding a constant to the entire magnitude spectrum.

Thank you vey much for your help!

Maria-Daphne
Hello Michael,

Thank you very much for your help! Below are my comments:

>What is your application? Are you trying to model this signal using
ARIMA-type >models?

My application is looking at geophone recordings of ambient noise
vibrations. In seismology, some theory for recovering the earths structure
relies on cross-correlating these ambient noise recordings. This theory
relies on the assumption that these noise recordings are for earthquakes
that are random in space and time, ie that the time series should be
stationary. This is not the case in what we record however. Moreover, the
recorded time series often has peaks at specific frequencies which
overwhelm the rest of the frequencies. What is then common to seismologists
is to take this signal prior to doing the cross-correlations and applying
some type of spectral normalization. This will in turn broaden the band of
noise signal in the cross-correlation and combat degradation caused by
monochromatic persistent sources.

>If you normalize in the manner you suggest, on a frequency-by-frequncy
basis, you >indeed remove the relative frequency relationships. This does
whiten the sequence, >but any intermediate information is lost. Could you
not simply simulate a white->noise signal directly?

This spectral normalization is commonly accomplished within the
seismological community I mean by prewhitening the signal. I have in this
context seen two prewhitening schemes: one is to add a constant over the
entire spectrum of the signal (similar to what Jeff suggested (I suggest
that the objective when whitening a signal would be to add energy at
frequencies *other* than those carrying
the signal's basic content, and such process would not result in a perfectly
flat magnitude response.) and the other is to take the magnitude spectrum,
and divide it by itself (to get one). However, most researchers never do
this substitution of the spectrum by a constant, but instead divide it by a
smoothed magnitude (what I have tried to do in my code). I suspect they
dont substitute by a constant because that would be changing the data too
much, ie one would loose too much information by changing completely the
spectrum instead of quieting the undesired overshoots in a frequency
spectrum. Or perhaps they dont do it because of another reason (again I
suspect?): given that the signal is bandlimited if you substitute its
spectrum with a constant that is like having a boxcar in the frequency
domain, so the signal perhaps would be smeared in time??

Which brings me to another question: as a student I remember being advised
to window my data (by blackmann hanning etc) prior to fft so to increase
spectral resolution (which I see as: if you dont window in the time domain,
you are multiplying by a boxcar, which is convolving with a sinc in the
frequency domain). However, why wasnt advised to window (by blackmann
hanning etc) in the frequency domain prior to ifft to improve time
resolution? Encorporated in this problem of prewhitening, shouldnt I be
dividing by the smoothed spectrum, then multiplying the magnitude with a
Blackman and then ifft? And leave the phase intact the entire time right?

As far as ARMA modeling I have seen in seismology applications in which
these models are used to simulate earthquake ground motions, however not for
the purpose of making a geophone record white. I wonder why they are not
used in that context I will read further into the literature you proposed.
Thank you again for your help! Maria-Daphne

Hello Jeff,

Thank you very much for your help. Below are some of my comments:

>If you make the magnitude perfectly flat, then after inverse FFT you would
have an >infinite impulse (Dirac Delta
>function) in the time domain.

That wouldnt be the case however, if you had a random phase right? So if
my original signal is random, and I only make the magnitude spectrum 1, but
dont change the phase, I shouldnt be getting an impulse correct?
>I suggest that the objective when whitening a signal would be to add energy
at >frequencies *other* than those carrying
>the signal's basic content, and such process would not result in a
perfectly flat >magnitude response.

Yes, I think this is right, that would be equivalent to adding a constant to
the entire magnitude spectrum.

Thank you vey much for your help!

Maria-Daphne