Reply by Maria-Daphne Mangriotis●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
Reply by mdma...@gmail.com●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
Reply by Jeff Brower●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
Reply by Jeff Brower●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
Reply by mdma...@gmail.com●October 6, 20102010-10-06
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).