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

# prewhitening a signal in matlab

Started by ●October 6, 2010

Reply by ●October 7, 20102010-10-07

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

> 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

Reply by ●October 7, 20102010-10-07

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 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

Reply by ●October 7, 20102010-10-07

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

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

Reply by ●October 7, 20102010-10-07

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

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