Discrete Hilbert Transform Bug?
Greg Heath
1. x(t) is real
2. n = length(x) is even
3. X = fft(x)
4. abs(X(n/2+1)) > 0
then
z = hilbert(x) IS NOT analytic.
The proof follows from
1. Zc = fftshift(fft(z))
2. Zc(1) = X(n/2+1)
Therefore the negative frequency spectrum does no vanish.
Comments are welcome.
Hope this helps.
Greg
robert bristow-johnson
he...@ll.mit.edu wrote on 09/17/2002 02:07:
> I posted this on comp.soft-sys.matlab and got no response. Perhaps I
> will get a better response here.
>
> Note that in MATLAB:
>
> 1. Indices go from 1 to N.
because that is the way God meant it to be.
and God said "Let there be
N
X[k] = SUM{ x[n] * exp(-j*2*pi*(n-1)*(k-1)/N) }
n=1
" and there was the frequency domain.
Oppenheim and Schafer and the DSP textbook authors are heretics to
contradict that by saying:
N-1
X[k] = SUM{ x[n] * exp(-j*2*pi*n*k/N) }
n=0
> 2. When N is even and fftshift is used, the components from N/2 + 1 to
> N are shifted to negative frequencies.
no. that is the case when fftshift() is *not* used (on the frequency domain
data). whether or not fftshift() is used on the time domain data, the
elements from X(N/2+1) to X(N) are the negative frequencies reflected up to
the upper half as per the periodic nature of the DFT.
> 3. MATLAB's hilbert(x) yields x+j*y where y is the classical hilbert
> transform of x.
i think that is true.
> he...@ll.mit.edu (Greg Heath) wrote in message
> news:<9e910d7a.02091...@posting.google.com>...
>> If
>>
>> 1. x(t) is real
>> 2. n = length(x) is even
>> 3. X = fft(x)
>> 4. abs(X(n/2+1)) > 0
>>
>> then
>>
>> z = hilbert(x) IS NOT analytic.
>>
>> The proof follows from
>>
>> 1. Zc = fftshift(fft(z))
>> 2. Zc(1) = X(n/2+1)
>>
>> Therefore the negative frequency spectrum does not vanish.
>>
>> Comments are welcome.
dunno.
--
r b-j
Wave Mechanics, Inc.
45 Kilburn St.
Burlington VT 05401-4750
tel: 802/951-9700 ext. 207 http://www.wavemechanics.com/
fax: 802/951-9799 rob...@wavemechanics.com
--
Greg Heath
> "Greg Heath" <he...@ll.mit.edu> schrieb im Newsbeitrag
> news:9e910d7a.02091...@posting.google.com...
> > I posted this on comp.soft-sys.matlab and got no response. Perhaps I
> > will get a better response here.
> >
> > Note that in MATLAB:
> >
> > 1. Indices go from 1 to N.
> > N are shifted to negative frequencies.
> > transform of x.
> >
> > news:<9e910d7a.02091...@posting.google.com>...
> > >
> > > 1. x(t) is real
> > > 2. n = length(x) is even
> > > 3. X = fft(x)
> > > 4. abs(X(n/2+1)) > 0
> > >
> > > then
> > >
> > > z = hilbert(x) IS NOT analytic.
Wrong.
z *IS* a discrete-time analytic signal.
See the references (1. O&S , 2. Marple) in
help hilbert
My confusion results fron associating the Nyquist component with the
negative freqency interval since the matlab function fftshift moves it
there.
> Vectors aren't analytic.
My guess is that you object to "analytic" instead of "analytic
signal".
If not, please explain.
> > > The proof follows from
> > >
> > > 1. Zc = fftshift(fft(z))
> > > 2. Zc(1) = X(n/2+1)
> > >
> > > Therefore the negative frequency spectrum does not vanish.
>
> Matlab is calculating the "analytical signal" of an input signal via
> discrete convolution rather than setting the negative frequency components
> to zero.
I never said Matlab *set* anything to zero. However, you would get the
same result if they did. See the algorithms in Marple.
Are you really trying to make a distinction between
X(N/2 + 2:N) = 0*X(N/2+2:N);
and
X(N/2 + 2:N) = 0;
?
Perplexed,
Greg
Greg Heath
> In article 9e910d7a.02091...@posting.google.com, Greg Heath at
> he...@ll.mit.edu wrote on 09/17/2002 02:07:
>
> > I posted this on comp.soft-sys.matlab and got no response. Perhaps I
> > will get a better response here.
> >
> > Note that in MATLAB:
> >
> > 1. Indices go from 1 to N.
>
> because that is the way God meant it to be.
>
> and God said "Let there be
>
> N
> X[k] = SUM{ x[n] * exp(-j*2*pi*(n-1)*(k-1)/N) }
> n=1
>
> " and there was the frequency domain.
>
> Oppenheim and Schafer and the DSP textbook authors are heretics to
> contradict that by saying:
>
> N-1
> X[k] = SUM{ x[n] * exp(-j*2*pi*n*k/N) }
> n=0
This is irrelevant w.r.t. the original post.
> > 2. When N is even and fftshift is used, the components from N/2 + 1 to
> > N are shifted to negative frequencies.
>
> no. that is the case when fftshift() is *not* used (on the frequency domain
> data). whether or not fftshift() is used on the time domain data, the
> elements from X(N/2+1) to X(N) are the negative frequencies reflected up to
> the upper half as per the periodic nature of the DFT.
I guess I didn't make my self clear:
When N is even, fft yields components with indices
[1,...,N/2+1,N/2+2,..., N] corresponding to freqencies
[0,..., N/2,N/2+1,...,N-1]*df
When fftshift is applied to the resulting fft, indices
[ 1,...,N/2,N/2+1,N/2+2,..., N] now correspond to previous
components
[ N/2+1,..., N, 1, 2,...,N/2] corresponding to freqencies
[-(N/2+1),..., -1, 0, 1,...,N/2]*df
For that reason, I associated the N/2+1 (Nyquist) component with
the negative frequency spectrum.
However, the hilbert transform references (O&S, Marple)associate the
Nyquist component with the positive frequency spectrum.
Hence, my confusion.
robert bristow-johnson
he...@ll.mit.edu wrote on 09/24/2002 13:41:
> robert bristow-johnson <rob...@wavemechanics.com> wrote in message
> news:<B9AF966D.77FC%rob...@wavemechanics.com>...
>> In article 9e910d7a.02091...@posting.google.com, Greg Heath at
>> he...@ll.mit.edu wrote on 09/17/2002 02:07:
>>
>>> I posted this on comp.soft-sys.matlab and got no response. Perhaps I
>>> will get a better response here.
>>>
>>> Note that in MATLAB:
>>>
>>> 1. Indices go from 1 to N.
>>
>> because that is the way God meant it to be.
>>
>> and God said "Let there be
>>
>> N
>> X[k] = SUM{ x[n] * exp(-j*2*pi*(n-1)*(k-1)/N) }
>> n=1
>>
>> " and there was the frequency domain.
>>
>> Oppenheim and Schafer and the DSP textbook authors are heretics to
>> contradict that by saying:
>>
>> N-1
>> X[k] = SUM{ x[n] * exp(-j*2*pi*n*k/N) }
>> n=0
>
> This is irrelevant w.r.t. the original post.
you think?
>>> 2. When N is even and fftshift is used, the components from N/2 + 1 to
>>> N are shifted to negative frequencies.
>>
>> no. that is the case when fftshift() is *not* used (on the frequency domain
>> data). whether or not fftshift() is used on the time domain data, the
>> elements from X(N/2+1) to X(N) are the negative frequencies reflected up to
>> the upper half as per the periodic nature of the DFT.
>
> I guess I didn't make my self clear:
>
> When N is even, fft yields components with indices
> [1,...,N/2+1,N/2+2,..., N] corresponding to freqencies
> [0,..., N/2,N/2+1,...,N-1]*df
(lessee: "1" is 0df and "N/2+1" is N/2*df, etc. and you were saying *what*
was irrelevant?) :-)
just FYI, Greg: it totally violates my sense of order when a very expensive
product that promises to use "familiar mathematical notation" and to express
"problems and solutions ... just as they are written mathematically" says,
OTOH that "1" is 0df and "N/2+1" is N/2*df and "N" is (N-1)*df. do you see
how you have to qualify your notation just to get the point across? if
Matlab lived up to it's promise, you would not have to do that, especially,
for such a simple and basic issue. Cleve and the MathWorks folk don't see
this as a problem, so each and every time i see someone getting messed up
with the notation simply because Matlab doesn't deal with indices that are
non-positive, i made it my mission in life to point it out to Cleve. that
the obstinately refuse to fix this (especially when i have spelled out how,
or what needs to get fixed), is inexcusable for any tool that has attained
such a status as a "standard".
> When fftshift is applied to the resulting fft,
okay, fine. we usually don't do that. we might apply fftshift to the
result of the ifft (and often to the input of the fft), but i haven't yet
seen it useful to apply fftshift to the result of the fft (i guess if you
want those negative frequencies to lie to the left of the positive one,
applying fftshift will do that. and since MATLAB messes up your DFT indices
(by 1) anyway, why not? it's just another housekeeping detail that you have
to remember.
> indices
> [ 1,...,N/2,N/2+1,N/2+2,..., N] now correspond to previous
> components
> [ N/2+1,..., N, 1, 2,...,N/2] corresponding to freqencies
> [-(N/2+1),..., -1, 0, 1,...,N/2]*df
and here, my "irrelevant" point gotcha! i think you meant to say:
When fftshift is applied to the resulting fft, indices
[ 1,...,N/2,N/2+1,N/2+2,..., N] now correspond to previous components
[ N/2+1,..., N, 1, 2,...,N/2 ] corresponding to frequencies
[-(N/2),..., -1, 0, 1,...,N/2-1]*df
did you not?
> For that reason, I associated the N/2+1 (Nyquist) component with
> the negative frequency spectrum.
sure, but it could just as well be positive. i look at it as a circular
counting issue precisely like associating 0x8000 with -32768 instead of
+32768 (for 16 bit int).
> However, the hilbert transform references (O&S, Marple)associate the
> Nyquist component with the positive frequency spectrum.
actually, since the DFT or DFS is used, they treat all of the frequencies as
"positive" since the spectrums are periodic anyway.
> Hence, my confusion.
okay, here's how to lose to confusion: set the nyquist component to zero.
it really should be there anyway.
if you have some non-zero component resulting from FFTing some real data,
that means there is a *cosine* component at Nyquist and that is the only one
you can represent (you cannot represent the sine component at Nyquist).
but the Hilbert transform of each cosine component is a sine. so you're
screwed, but the Nyquist/Shannon Sampling Theorem predicted that anyway.
you must sample at a rate that exceeds twice your highest frequency
component.
like DC, the Hilbert transform will map any Nyquist component to zero.
Greg Heath
> In article 9e910d7a.02092...@posting.google.com, Greg Heath at
> he...@ll.mit.edu wrote on 09/24/2002 13:41:
>
> > robert bristow-johnson <rob...@wavemechanics.com> wrote in message
> > news:<B9AF966D.77FC%rob...@wavemechanics.com>...
> >> In article 9e910d7a.02091...@posting.google.com, Greg Heath at
> >> he...@ll.mit.edu wrote on 09/17/2002 02:07:
> >>
> >>> I posted this on comp.soft-sys.matlab and got no response. Perhaps I
> >>> will get a better response here.
> >>>
> >>> Note that in MATLAB:
> >>>
> >>> 1. Indices go from 1 to N.
> >>
> >> because that is the way God meant it to be.
> >>
> >> and God said "Let there be
> >>
> >> N
> >> X[k] = SUM{ x[n] * exp(-j*2*pi*(n-1)*(k-1)/N) }
> >> n=1
> >>
> >> " and there was the frequency domain.
> >>
> >> Oppenheim and Schafer and the DSP textbook authors are heretics to
> >> contradict that by saying:
> >>
> >> N-1
> >> X[k] = SUM{ x[n] * exp(-j*2*pi*n*k/N) }
> >> n=0
> >
> > This is irrelevant w.r.t. the original post.
>
> you think?
Sure. Regardless of which convention is used, is the Nyquist
component
considered to be in the positive frequency range or the negative
frequency range? The hilbert transform references O&S and Marple ==>
+,
while I, using fftshift,was led to -.
> >>> 2. When N is even and fftshift is used, the components from N/2 + 1 to
> >>> N are shifted to negative frequencies.
> >>
> >> no. that is the case when fftshift() is *not* used (on the frequency domain
> >> data). whether or not fftshift() is used on the time domain data, the
> >> elements from X(N/2+1) to X(N) are the negative frequencies reflected up to
> >> the upper half as per the periodic nature of the DFT.
> >
> > I guess I didn't make my self clear:
> >
> > When N is even, fft yields components with indices
> > [1,...,N/2+1,N/2+2,..., N] corresponding to freqencies
> > [0,..., N/2,N/2+1,...,N-1]*df
>
> (lessee: "1" is 0df and "N/2+1" is N/2*df, etc. and you were saying *what*
> was irrelevant?) :-)
Irrelevant to solving my immediate problem.
> just FYI, Greg: it totally violates my sense of order when a very expensive
> product that promises to use "familiar mathematical notation" and to express
> "problems and solutions ... just as they are written mathematically" says,
> OTOH that "1" is 0df and "N/2+1" is N/2*df and "N" is (N-1)*df. do you see
> how you have to qualify your notation just to get the point across? if
> Matlab lived up to it's promise, you would not have to do that, especially,
> for such a simple and basic issue. Cleve and the MathWorks folk don't see
> this as a problem, so each and every time i see someone getting messed up
> with the notation simply because Matlab doesn't deal with indices that are
> non-positive, i made it my mission in life to point it out to Cleve. that
> the obstinately refuse to fix this (especially when i have spelled out how,
> or what needs to get fixed), is inexcusable for any tool that has attained
> such a status as a "standard".
I understand your frustration. However I was weaned on Fortran in the
early 60's so I've long since adjusted to that. My understanding is
that originally MATLAB started with a Fortran underpinning and got
locked into the 1 to N indexing. However, since later Fortran versions
accommodate more flexible indexing, I agree that MATLAB should also.
> > When fftshift is applied to the resulting fft,
>
> okay, fine. we usually don't do that. we might apply fftshift to the
> result of the ifft (and often to the input of the fft), but i haven't yet
> seen it useful to apply fftshift to the result of the fft
???
But that seems to be exactly what is recommended in the inept
discussions of 'help fftshift' and 'help ifftshift'! They imply the
proper usage is
X = fft(x);
Xc = fftshift(X); % to "center" the DC frequency
X = ifftshift(Xc); % to "undo" fftshift
Anyway, when N is even, fftshift and ifftshift are equivalent, so that
is not the answer to my immediate problem. For example,
X = [1 2 3 4 5 6 7 8] % The unique components are 1-5
fftshift(X) = [5 6 7 8 1 2 3 4] % The Nyquist component is "5"
ifftshift(X) = [5 6 7 8 1 2 3 4] % The Nyquist component is mapped
% to the negative
frequency interval
Although the original post is only concerned with N even, notice that,
when N is odd, the shift results are different:
X = [1 2 3 4 5 6 7] % The unique components are 1-4
Xc = fftshift(X) = [5 6 7 1 2 3 4] % The Nyquist component is "4.5"
X = ifftshift(Xc) = [1 2 3 4 5 6 7]
whereas
Y = ifftshift(X) = [4 5 6 7 1 2 3]
X = fftshift(Y) = [1 2 3 4 5 6 7]
So I assume you always use N even.
> ... (i guess if you
> want those negative frequencies to lie to the left of the positive one,
> applying fftshift will do that. and since MATLAB messes up your DFT indices
> (by 1) anyway, why not? it's just another housekeeping detail that you have
> to remember.
>
> > indices
> > [ 1,...,N/2,N/2+1,N/2+2,..., N] now correspond to previous
> > components
> > [ N/2+1,..., N, 1, 2,...,N/2] corresponding to freqencies
> > [-(N/2+1),..., -1, 0, 1,...,N/2]*df
>
> and here, my "irrelevant" point gotcha! i think you meant to say:
>
> When fftshift is applied to the resulting fft, indices
> [ 1,...,N/2,N/2+1,N/2+2,..., N] now correspond to previous components
> [ N/2+1,..., N, 1, 2,...,N/2 ] corresponding to frequencies
> [-(N/2),..., -1, 0, 1,...,N/2-1]*df
>
> did you not?
No. I didn't.
> > For that reason, I associated the N/2+1 (Nyquist) component with
> > the negative frequency spectrum.
>
> sure, but it could just as well be positive. i look at it as a circular
> counting issue precisely like associating 0x8000 with -32768 instead of
> +32768 (for 16 bit int).
I agree with your point of view. However, the question at hand is to
decipher the apparent discrepancy between MATLAB and it's references
(O&S, Marple):
Is the Nyquist component considered to be in the positive or negative
part of the spectrum?
> > However, the hilbert transform references (O&S, Marple)associate the
> > Nyquist component with the positive frequency spectrum.
>
> actually, since the DFT or DFS is used, they treat all of the frequencies as
> "positive" since the spectrums are periodic anyway.
No!
Please refer to Marple's paper and Ch. 10 in O&S w.r.t. the hilbert
transform and positive/negative parts of the spectrum.
> > Hence, my confusion.
>
> okay, here's how to lose to confusion: set the nyquist component to zero.
> it really should be there anyway.
Agree if you mean "it really should *not* be there anyway."
> if you have some non-zero component resulting from FFTing some real data,
> that means there is a *cosine* component at Nyquist and that is the only one
> you can represent (you cannot represent the sine component at Nyquist).
Agree.
> but the Hilbert transform of each cosine component is a sine. so you're
> screwed, but the Nyquist/Shannon Sampling Theorem predicted that anyway.
> you must sample at a rate that exceeds twice your highest frequency
> component.
>
> like DC, the Hilbert transform will map any Nyquist component to zero.
Agree.
But that still leaves the inconsistency that O&S and Marple consider
the Nyquist component to be in the + part of the spectrum whereas
MATLAB appears to consider it in the - part.
My question:
Is one of these interpretations considered "right" by convention?
or
Both are acceptable and you have to choose depending on what you are
doing.
By the way, I'm trying to use the hilbert transform to obtain a causal
time domain response from (only) the sampled amplitude of the
continuous-time Fourier Transform. When N is odd things work o.k.
When N is even, the Nyquist component yields spurious oscillations in
the hilbert transform phase estimate.
Currently, I'm just trying to understand, in more detail, what's going
on.
Thanks for your help.
Greg
robert bristow-johnson
Greg,
i meant to respond to you directly when i got your email but did not get
around to it. my apologies.
In article 9e910d7a.0209...@posting.google.com, Greg Heath at
the Nyquist component is considered to be both the + and -. for a real
signal the spectrum must be symmetric about 0 (even symmetry for real part
and magnitude, odd symmetry for imaginary part and phase). in order for a
non-zero Nyquist component to result in a real sinusoid in the
continuous-time domain, 1/2 of its amplitude is placed at -Nyquist and 1/2
at +Nyquist. and the result for the Nyquist component will be a pure
cosine, no sine. when you sample back to the discrete world, you get the
original discrete samples and discrete DFT. in the DFT frequency domain,
the two 1/2 amplitude components (which, like DC, must be real if the
original signal was ever real) will overlap-add (because of imagining) be in
phase and add to the single Nyquist amplitude. if there ever was a sine
component at Nyquist, any information regarding that component is lost when
it was sampled or hypothetically sampled.
anyway, this is one reason i do not like any non-zero amplitude in the
Nyquist bin. i think *all* of my FFT based analysis/synthesis code
deliberately sets both the Nyquist bin and the DC bin to zero just so i
don't have to deal with anything weird.
>>>>> 2. When N is even and fftshift is used, the components from N/2 + 1 to
>>>>> N are shifted to negative frequencies.
>>>>
>>>> no. that is the case when fftshift() is *not* used (on the frequency
>>>> domain
>>>> data). whether or not fftshift() is used on the time domain data, the
>>>> elements from X(N/2+1) to X(N) are the negative frequencies reflected up to
>>>> the upper half as per the periodic nature of the DFT.
>>>
>>> I guess I didn't make my self clear:
>>>
>>> When N is even, fft yields components with indices
>>> [1,...,N/2+1,N/2+2,..., N] corresponding to freqencies
>>> [0,..., N/2,N/2+1,...,N-1]*df
>>
>> (lessee: "1" is 0df and "N/2+1" is N/2*df, etc. and you were saying *what*
>> was irrelevant?) :-)
>
> Irrelevant to solving my immediate problem.
perhaps irrelevant, but having to make that point, that "1" is 0*df, "2" is
1*df, ... "N/2+1" is N/2*df is an unnecessary pain in the rear that makes my
point regarding MATLAB.
wow! i didn't even know there *was* an "ifftshift()". for 1-dim vectors,
fftshift() undoes itself. for this context, fftshift() simply swaps the 1st
half and the 2nd half.
> Anyway, when N is even, fftshift and ifftshift are equivalent,
i never thought about the N odd case since i always did FFTs with N = 2^p.
i tried a simple fftshift on an odd N and
誌=[1 2 3 4 5 6 7];
蜻ftshift(x)
ans =
5 6 7 1 2 3 4
so it looks like MATLAB sees the just below Nyquist bin (#4) as positive as
it should and put the DC bin (#1) in the middle.
> so that
> is not the answer to my immediate problem. For example,
>
> X = [1 2 3 4 5 6 7 8] % The unique components are 1-5
> fftshift(X) = [5 6 7 8 1 2 3 4] % The Nyquist component is "5"
> ifftshift(X) = [5 6 7 8 1 2 3 4] % The Nyquist component is mapped
> % to the negative
> frequency interval
>
> Although the original post is only concerned with N even, notice that,
> when N is odd, the shift results are different:
>
> X = [1 2 3 4 5 6 7] % The unique components are 1-4
> Xc = fftshift(X) = [5 6 7 1 2 3 4] % The Nyquist component is "4.5"
> X = ifftshift(Xc) = [1 2 3 4 5 6 7]
>
> whereas
>
> Y = ifftshift(X) = [4 5 6 7 1 2 3]
> X = fftshift(Y) = [1 2 3 4 5 6 7]
>
> So I assume you always use N even.
yes i have. in fact i always used N = power of 2.
anyway, i use fftshift() to fix a segment of windowed audio (so the bulge is
in the middle, around sample N/2) so that the center of the window (and the
bulge of the windowed audio) is centered at t=0. then, when i FFT this, the
FFT looks at it sorta like i have an undelayed segment of audio centered at
zero. if i don't do this, then the phase of every other bin in the FFT
result has this messy jump of pi or -pi radians. if i *do* this preswap of
time domain data before FFT gets it, then my phase response is smooth. then
i mess around with the spectrum and inverse FFT back to time and then i have
to fix my original fftshift() with another one that undoes it. i will have
to look into using ifftshift() to make the code more bullet proof in case i
ever use an odd N.
>> ... (i guess if you
>> want those negative frequencies to lie to the left of the positive one,
>> applying fftshift will do that. and since MATLAB messes up your DFT indices
>> (by 1) anyway, why not? it's just another housekeeping detail that you have
>> to remember.
>>
>>> indices
>>> [ 1,...,N/2,N/2+1,N/2+2,..., N] now correspond to previous components
>>> [ N/2+1,..., N, 1, 2,...,N/2] corresponding to freqencies
>>> [-(N/2+1),..., -1, 0, 1,...,N/2]*df
>>
>> and here, my "irrelevant" point gotcha! i think you meant to say:
>>
>> When fftshift is applied to the resulting fft, indices
>> [ 1,...,N/2,N/2+1,N/2+2,..., N] now correspond to previous components
>> [ N/2+1,..., N, 1, 2,...,N/2 ] corresponding to frequencies
>> [-(N/2),..., -1, 0, 1,...,N/2-1]*df
>>
>> did you not?
>
> No. I didn't.
of course you did, Greg. if you meant that bin "1" corresponds to frequency
0*df and bin "2" corresponds to frequency 1*df, how could you consistently
intend to say bin "N/2" corresponds to frequency N/2*df or bin "N/2+1"
corresponds to frequency -(N/2+1)*df? of course, in MATLAB, bin "N/2"
corresponds to frequency (N/2-1)*df and bin "N/2+1" corresponds to frequency
(N/2)*df or, if you fold it back frequency -(N/2)*df.
you are underscoring my point that MATLAB insistence to assign frequencies
to bins that are misnumbered by 1 does not serve as a conduit to clear and
error free discussion of the arcane technical issues. this old mistake of
MATLAB's gets in the way and increases confusion. that is why it is simply
wrong and should get fixed.
>>> For that reason, I associated the N/2+1 (Nyquist) component with
>>> the negative frequency spectrum.
>>
>> sure, but it could just as well be positive. i look at it as a circular
>> counting issue precisely like associating 0x8000 with -32768 instead of
>> +32768 (for 16 bit int).
>
> I agree with your point of view. However, the question at hand is to
> decipher the apparent discrepancy between MATLAB and it's references
> (O&S, Marple):
>
> Is the Nyquist component considered to be in the positive or negative
> part of the spectrum?
assuming that the discrete samples that went into the FFT were sampled from
a continuous-time signal that conceptually had frequencies not limited to
any bandwidth, then half of the Nyquist component is positive and half is
negative.
>>> However, the hilbert transform references (O&S, Marple)associate the
>>> Nyquist component with the positive frequency spectrum.
>>
>> actually, since the DFT or DFS is used, they treat all of the frequencies as
>> "positive" since the spectrums are periodic anyway.
>
> No!
>
> Please refer to Marple's paper and Ch. 10 in O&S w.r.t. the hilbert
> transform and positive/negative parts of the spectrum.
i don't have Marple but i can't agree with you about O&S. say Eq. (10.39).
those are positive n that get set the x[n] to zero. but because of
periodicity (with N being the period), "N/2 < n < N" is essentially the same
as saying "-N/2 < n < 0", anyway.
>>> Hence, my confusion.
>>
>> okay, here's how to lose to confusion: set the nyquist component to zero.
>> it really should be there anyway.
>
> Agree if you mean "it really should *not* be there anyway."
yup. (oops ... typo)
>> if you have some non-zero component resulting from FFTing some real data,
>> that means there is a *cosine* component at Nyquist and that is the only one
>> you can represent (you cannot represent the sine component at Nyquist).
>
> Agree.
>
>> but the Hilbert transform of each cosine component is a sine. so you're
>> screwed, but the Nyquist/Shannon Sampling Theorem predicted that anyway.
>> you must sample at a rate that exceeds twice your highest frequency
>> component.
>>
>> like DC, the Hilbert transform will map any Nyquist component to zero.
>
> Agree.
>
> But that still leaves the inconsistency that O&S and Marple consider
> the Nyquist component to be in the + part of the spectrum whereas
> MATLAB appears to consider it in the - part.
>
> My question:
>
> Is one of these interpretations considered "right" by convention?
no. neither can be more "right" than the other.
> or
>
> Both are acceptable and you have to choose depending on what you are
> doing.
if there is *any* non-zero amplitude in the Nyquist bin, you have to split
it half and half to correspond to a *real* signal. but since you're doing
Hilbert stuff anyway, my suspicion is that you are messing around with
"analytical signals" which are complex. in that case, then the Hilbert
component will be positive (all of the negative frequencies are set to zero
amplitude), BUT you will have a problem. the real part of your analytic
signal will have this non-zero cosine and the imag part will have no
corresponding sine. so, again, i would nuke the Nyquist bin. explicitly
annihilate! set that sucker to zero.
> By the way, I'm trying to use the hilbert transform to obtain a causal
> time domain response from (only) the sampled amplitude of the
> continuous-time Fourier Transform. When N is odd things work o.k.
> When N is even, the Nyquist component yields spurious oscillations in
> the hilbert transform phase estimate.
well, you still have assumptions to make. will it be minimum-phase or
something else. there are many causal impulse responses that correspond to
any single magnitude frequency response.
anyway you can take your magnitude response, attach any phase you like,
inverse FFT and use fftshift() on that and you have a causal impulse
response. if you would like the last half of the impulse response to be
zero, then the real and imag parts of the frequency response have to be a
Hilbert transform pair.
> Currently, I'm just trying to understand, in more detail, what's going
> on.
>
> Thanks for your help.
FWIW, sorry for the delay.
--
robert bristow-johnson
rob...@wavemechanics.com wrote on 09/26/2002 20:33:
> In article 9e910d7a.0209...@posting.google.com, Greg Heath at
> he...@ll.mit.edu wrote on 09/26/2002 15:39:
...
>
>> By the way, I'm trying to use the hilbert transform to obtain a causal
>> time domain response from (only) the sampled amplitude of the
>> continuous-time Fourier Transform. When N is odd things work o.k.
>> When N is even, the Nyquist component yields spurious oscillations in
>> the hilbert transform phase estimate.
>
> well, you still have assumptions to make. will it be minimum-phase or
> something else. there are many causal impulse responses that correspond to
> any single magnitude frequency response.
>
> anyway you can take your magnitude response, attach any phase you like,
> inverse FFT and use fftshift() on that and you have a causal impulse
> response. if you would like the last half of the impulse response to be
> zero, then the real and imag parts of the frequency response have to be a
> Hilbert transform pair.
BTW, i meant to add that if your problem boils down to, you are given a
magnitude frequency response:
|H(w)|
where
H(w) = Hr(w) + j*Hi(w)
and you want a solution (for Hr and Hi) to:
Hr(w)^2 + Hi(w)^2 = |H(w)|^2 (the square of the given)
and
Hi(w) = -Hilbert{ Hr(w) } (which makes h(t) causal)
then i don't have an answer, but i am thinking about it with idle CPU
cycles.
i dunno if that is the issue or not.
L8r,
r b-j
Paavo Jumppanen
.snip
.
>
> > but the Hilbert transform of each cosine component is a sine. so you're
> > screwed, but the Nyquist/Shannon Sampling Theorem predicted that anyway.
> > you must sample at a rate that exceeds twice your highest frequency
> > component.
> >
> > like DC, the Hilbert transform will map any Nyquist component to zero.
>
> Agree.
>
> But that still leaves the inconsistency that O&S and Marple consider
> the Nyquist component to be in the + part of the spectrum whereas
> MATLAB appears to consider it in the - part.
>
> My question:
>
> Is one of these interpretations considered "right" by convention?
>
> or
>
> Both are acceptable and you have to choose depending on what you are
> doing.
In the digital domain I would argue that both could be considered
valid. Remember that a sampled time signal has a periodic spectrum. In
the z-domain the Nyquist sample is purely real and lies immediately
between the +ve and -ve spectrum parts (ie. unit circle in z-domain)
as with the dc component. This is in contrast with the continuous time
case where the +ve and -ve spectrum parts only ever meet at DC.
Regards,
Paavo Jumppanen
Author of AtSpec : A 2 channel PC based FFT Spectrum Analyzer
http://www.taquis.com
Greg Heath
> Date: Thu, 26 Sep 2002 20:33:47 -0400
> From: robert bristow-johnson <rob...@wavemechanics.com>
> To: he...@ll.mit.edu
> Cc: Greg Heath <he...@ll.mit.edu>
> Newgroups: comp.dsp, comp.soft-sys.matlab
> Subject: Re: Discrete Hilbert Transform Bug? (Cleve, look at this.)
>
> Greg,
>
> i meant to respond to you directly when i got your email but did not
> get around to it. my apologies.
None needed. But I'll accept it anyway (I get so few of them).
Of course. Just like the DC component.
> for a real
> signal the spectrum must be symmetric about 0 (even symmetry for real part
> and magnitude, odd symmetry for imaginary part and phase). in order for a
> non-zero Nyquist component to result in a real sinusoid in the
> continuous-time domain, 1/2 of its amplitude is placed at -Nyquist and 1/2
> at +Nyquist. and the result for the Nyquist component will be a pure
> cosine, no sine. when you sample back to the discrete world, you get the
> original discrete samples and discrete DFT. in the DFT frequency domain,
> the two 1/2 amplitude components (which, like DC, must be real if the
> original signal was ever real) will overlap-add (because of imagining) be in
> phase and add to the single Nyquist amplitude. if there ever was a sine
> component at Nyquist, any information regarding that component is lost when
> it was sampled or hypothetically sampled.
Agree.
> anyway, this is one reason i do not like any non-zero amplitude in the
> Nyquist bin. i think *all* of my FFT based analysis/synthesis code
> deliberately sets both the Nyquist bin and the DC bin to zero just so i
> don't have to deal with anything weird.
I wanted to but was afraid that I would just add to my confusion. so...
I tried the Bilinear Transform before the Hilbert Transform. That works
for obtaining Hi(w) from Hr(w) when N is arbitrary or for obtaining
phase(w) from log(|H(w)|) when N is odd (when N is even, log-amp -> -oo
for the Nyqist component).
So, I decided that, If I had to, I would only use N odd. Then I would
investigate anti-aliasing frequency filtering in order to use N even.
> >>>>> 2. When N is even and fftshift is used, the components from N/2 + 1 to
> >>>>> N are shifted to negative frequencies.
> >>>>
> >>>> no. that is the case when fftshift() is *not* used (on the frequency
> >>>> domain
> >>>> data). whether or not fftshift() is used on the time domain data, the
> >>>> elements from X(N/2+1) to X(N) are the negative frequencies reflected up
> >>>> to the upper half as per the periodic nature of the DFT.
> >>>
> >>> I guess I didn't make my self clear:
> >>>
> >>> When N is even, fft yields components with indices
> >>> [1,...,N/2+1,N/2+2,..., N] corresponding to freqencies
> >>> [0,..., N/2,N/2+1,...,N-1]*df
> >>
> >> (lessee: "1" is 0df and "N/2+1" is N/2*df, etc. and you were saying *what*
> >> was irrelevant?) :-)
> >
> > Irrelevant to solving my immediate problem.
>
> perhaps irrelevant, but having to make that point, that "1" is 0*df, "2" is
> 1*df, ... "N/2+1" is N/2*df is an unnecessary pain in the rear that makes my
> point regarding MATLAB.
Agree.
Right. N odd works fine (except for the obvious loss in speed).
It's not my fault! The MATLAB indexing ...
> you are underscoring my point that MATLAB insistence to assign frequencies
> to bins that are misnumbered by 1 does not serve as a conduit to clear and
> error free discussion of the arcane technical issues. this old mistake of
> MATLAB's gets in the way and increases confusion. that is why it is simply
> wrong and should get fixed.
O.K. add another signature to the petition.
> >>> For that reason, I associated the N/2+1 (Nyquist) component with
> >>> the negative frequency spectrum.
> >>
> >> sure, but it could just as well be positive. i look at it as a circular
> >> counting issue precisely like associating 0x8000 with -32768 instead of
> >> +32768 (for 16 bit int).
> >
> > I agree with your point of view. However, the question at hand is to
> > decipher the apparent discrepancy between MATLAB and it's references
> > (O&S, Marple):
> >
> > Is the Nyquist component considered to be in the positive or negative
> > part of the spectrum?
>
> assuming that the discrete samples that went into the FFT were sampled from
> a continuous-time signal that conceptually had frequencies not limited to
> any bandwidth, then half of the Nyquist component is positive and half is
> negative.
I agree. But when you sample to get the discrete spectrum an alternate
approach is required.
> >>> However, the hilbert transform references (O&S, Marple)associate the
> >>> Nyquist component with the positive frequency spectrum.
> >>
> >> actually, since the DFT or DFS is used, they treat all of the frequencies
> >> as "positive" since the spectrums are periodic anyway.
> >
> > No!
> >
> > Please refer to Marple's paper and Ch. 10 in O&S w.r.t. the hilbert
> > transform and positive/negative parts of the spectrum.
>
> i don't have Marple but i can't agree with you about O&S. say Eq. (10.39).
> those are positive n that get set the x[n] to zero. but because of
> periodicity (with N being the period), "N/2 < n < N" is essentially the same
> as saying "-N/2 < n < 0", anyway.
Then the centered spectrum consists of 4 parts:
1. Nyquist frequency
2. negative frequencies
3. zero frequency
4. positive frequencies
Take the 'help hilbert' example with Fs = 4:
f = [ 0 1 2 3];
fc = [ -2 -1 0 1];
Question: neg? neg zero pos.
Answer: No
x = [ 1 2 3 4];
z = hilbert(x) = [1+1i 2-1i 3-1i 4];
Z = fft(z) = [ 10 -4+4i -2 0];
Zc = fftshift(Z) = [ -2 0 10 -4+4i];
Nyq. neg. zero pos.
> >>> Hence, my confusion.
> >>
> >> okay, here's how to lose to confusion: set the nyquist component to zero..
With N odd and assuming the transform is an analytic signal in frequency,
hilbert fits the bill. I'm dealing with troposheric and ionospheric
noninteger power law spectra. So the typical rational fraction DSP
techniques are not applicable.
> anyway you can take your magnitude response, attach any phase you like,
> inverse FFT and use fftshift() on that and you have a causal impulse
> response.
Not always. I have come up with U shaped impulse responses that don't be
come causal with fftshift.
> if you would like the last half of the impulse response to be
> zero, then the real and imag parts of the frequency response have to be a
> Hilbert transform pair.
Yeah. Unfortunately, that doesn't happen because my log-amplitude and
phase are the hilbert transform pair! Nevertheless, I do get causal
responses like t^v exp(-w0 t) u(t).
Greg