Analytic Signal Generation in the Frequency Domain

Hi All,

I've been trying to generate an analytic signal (Hilbert
transformation) in the frequency domain using the method outlined in
the IEEE Marple paper in the Transactions on Signal Processing
"Computing the Discrete-Time Analytic Signal via FFT" (9/1999). I'm
also using Rick's UDSP 2nd edition book as a reference, which outlines
the same method.

The problem I'm haivng is that the resulting imaginary time domain
sequence has what appear to be severe errors unless specific conditions
are met on the input. However, the real sequence is correct.

To be more specific, here is pseudocode of the method I am using (x
contains an infinitely long real input sequence which will be processed
in blocks of size N). Overlap-and-add is used because circular
convolution appears to occur for this operation:

1. zero-pad x to length 2N
2. X = fft(x)
3. Perform Hilbert Transformation calculations in frequency domain
(outlined in Marple paper and Rick's book):
X[DC] *= 1
X[positive freqs] *= 2
X[Nyquist] *= 1
X[negative freqs] = 0
4. y = ifft(X)
5. Overlap and add y_previous (y[N:2N-1] from previous block) to
y[0:N-1]
6. Store y[N:2N-1] to y_previous

The real part of the resulting y is always correct. However, the
imaginary part is only correct if the real input x is a sinusoid with
either an integer period within the block or an integer period plus
half of a cycle. For example, the following sinusoids work correctly
when a block size of 64 is used and i is the index:

x=sin(2*pi*(1/64)*i)
x=sin(2*pi*(1.5/64)*i)

whereas this results in significant errors (frequency of 1.25/64):

x=sin(2*pi*(1.25/64)*i)

So here are my questions:

1. Is the above method (generating the analytic signal in the frequency
domain) supposed to work in block processing? If so, is overlap-and-add
required?
2. Is this method supposed to work for any arbitrary sinusoidal input?
3. If yes to all of the above, based on the pseudocode is there
something I might be doing wrong?

I notice that if I use this method for a long input (say several
seconds of audio) and take an FFT of the entire sample, it seems to
work except for oddities at the ends. But I have read references that
say this should work for block-based processing. Thanks for any help.

Regards,
Bob

BobM wrote:
> Hi All,
>
> I've been trying to generate an analytic signal (Hilbert
> transformation) in the frequency domain using the method outlined in
> the IEEE Marple paper in the Transactions on Signal Processing
> "Computing the Discrete-Time Analytic Signal via FFT" (9/1999). I'm
> also using Rick's UDSP 2nd edition book as a reference, which outlines
> the same method.
>
> The problem I'm haivng is that the resulting imaginary time domain
> sequence has what appear to be severe errors unless specific conditions
> are met on the input. However, the real sequence is correct.
>
> To be more specific, here is pseudocode of the method I am using (x
> contains an infinitely long real input sequence which will be processed
> in blocks of size N). Overlap-and-add is used because circular
> convolution appears to occur for this operation:
>
> 1. zero-pad x to length 2N
> 2. X = fft(x)
> 3. Perform Hilbert Transformation calculations in frequency domain
> (outlined in Marple paper and Rick's book):
> X[DC] *= 1
> X[positive freqs] *= 2
> X[Nyquist] *= 1
> X[negative freqs] = 0

this is multiplying your frequency response by a certain function.
what is the impulse response that corresponds to it?

> 4. y = ifft(X)
> 5. Overlap and add y_previous (y[N:2N-1] from previous block) to
> y[0:N-1]
> 6. Store y[N:2N-1] to y_previous

Bob, can you get away with OLA with this?  do you know ifdoing this
Hilbert Transform the way you do it is the same as convolving with an
impulse response that is short enough to do OLA?  (of length N or
smaller, with at least N zeros appended?)

r b-j

robert bristow-johnson wrote:
> BobM wrote:
> > 1. zero-pad x to length 2N
> > 2. X = fft(x)
> > 3. Perform Hilbert Transformation calculations in frequency domain
> > (outlined in Marple paper and Rick's book):
> > X[DC] *= 1
> > X[positive freqs] *= 2
> > X[Nyquist] *= 1
> > X[negative freqs] = 0
>
> this is multiplying your frequency response by a certain function.
> what is the impulse response that corresponds to it?

That's a very good question, but to be honest I haven't considered the
time domain characteristics of this function.

> > 4. y = ifft(X)
> > 5. Overlap and add y_previous (y[N:2N-1] from previous block) to
> > y[0:N-1]
> > 6. Store y[N:2N-1] to y_previous
>
> Bob, can you get away with OLA with this?  do you know ifdoing this
> Hilbert Transform the way you do it is the same as convolving with an
> impulse response that is short enough to do OLA?  (of length N or
> smaller, with at least N zeros appended?)

Another good question! Well, it looked as if OLA would be needed,
because it fixes the imaginary output when the input is an integer plus
half-cycle period within the block. But it doesn't solve the problem
for other cases. But honestly, I do not know if this is the same as
convolving with an impulse response that is short enough for OLA.

I have some Octave code that shows the problem I'm seeing if that is
helpful.

Regards,
Bob

Bob, i won't have time to deal with your octave code.  but i do want to
say that the DFT (of which the FFT is a fast method of doing it)
fundamentally maps a periodic sequence of period N (so you only need N
adjacent samples to completely describe it) from some domain (call it
the "time domain") to another periodic sequence of the same period in
the reciprocal domain (call that one the "frequency domain") and the
inverse DFT maps it back.  in this sense, it is accurate (but sometimes
disputed, we get in fights about this language) to say that the DFT
periodically extends the finite set of data passed to it.  the DFT
"thinks" the N samples passed to or returned from it are N adjacent
samples of a period sequence of period N even if those N samples came
from something else.

when multiplication is happening in one of the domains (that is what
you are doing when you compute the Analytic Signal the way you do) that
is equivalent to performing *circular* convolution in the other domain.
 (this "circular convolution" can be thought of as linear convolution
applied to the periodically extended data passed to the DFT.  so the
x[N-1] sample gets inherently followed by x[0].  in fact, *every*
x[N-1] sample gets followed immediately by x[0] and every x[0] is
preceded by x[N-1] forever in both left and right directions.

so the effect is that, if your DFT is of length N (you called it "2N",
but i don't want to use that convention because i'll get something
wrong), then if you had a zero padded piece of input data x[n] that was
of non-zero length of M (that is x[M], x[M+1], ... x[N-1] were all
zero) and you convolved it with an "impulse response" h[n] of length L
(so h[L], h[L+1], ... h[N-1] are all zero), then, if M and L are small
enough, the result of the convolution would have non-zero length of
M+L-1 and if that is less than or equal to N, the results of the
circular convolution and the linear convolution will agree for indices
0 to N-1.

if you're doing overlap-add, this means, if your impulse response is of
length L, the most samples you can get out of this using circular
convolution is N-L+1 for each overlapping frame. so you have to zero
pad more than just your input, you have to zero pad your impulse
response as well.

r b-j

On 15 Jan 2007 18:00:11 -0800, "BobM" <BobM.DSP@gmail.com> wrote:

>robert bristow-johnson wrote:
>> BobM wrote:
>> > 1. zero-pad x to length 2N
>> > 2. X = fft(x)
>> > 3. Perform Hilbert Transformation calculations in frequency domain
>> > (outlined in Marple paper and Rick's book):
>> > X[DC] *= 1
>> > X[positive freqs] *= 2
>> > X[Nyquist] *= 1
>> > X[negative freqs] = 0
>>
>> this is multiplying your frequency response by a certain function.
>> what is the impulse response that corresponds to it?

Hi Bob,

>That's a very good question, but to be honest I haven't considered the
>time domain characteristics of this function.

Yep, R B-J has a nasty habit of asking tough (but good) 
questions.  Because the freq-domain characteristic is 
rectangular, something like:

   1,2,2,2,...2,2,2,1,0,0,0,0,...,0,0,0,0

I'd guess the time-domain impulse response would 
some sort of complex sequence whose magnitude response 
would be very similar to a sinx(x) sequence.

>
>> > 4. y = ifft(X)
>> > 5. Overlap and add y_previous (y[N:2N-1] from previous block) to
>> > y[0:N-1]
>> > 6. Store y[N:2N-1] to y_previous
>>
>> Bob, can you get away with OLA with this?  do you know ifdoing this
>> Hilbert Transform the way you do it is the same as convolving with an
>> impulse response that is short enough to do OLA?  (of length N or
>> smaller, with at least N zeros appended?)
>
>Another good question! Well, it looked as if OLA would be needed,
>because it fixes the imaginary output when the input is an integer plus
>half-cycle period within the block. But it doesn't solve the problem
>for other cases. But honestly, I do not know if this is the same as
>convolving with an impulse response that is short enough for OLA.
>
>I have some Octave code that shows the problem I'm seeing if that is
>helpful.

I'm sorry to say, Bob, but I haven't modeled your 
process of OLA in conjunction with the "freq-domain 
analytic sequence generation' scheme.  
I'll do what I can to model (using MATLAB) this 
whole idea, and get back to you.  Perhaps, together, 
we'll figure this all out.

By the way Bob,
I'll be happy to send you an errata for
the 2nd edition of my book if you can tell 
me exactly what is the "Printing Number" of your 
copy of the book.  You can find your 
"Printing Number" on the page just before the 
"Dedication" page.

On that page (before the Dedication) you'll see 
all sorts of publisher-related information.
Down toward the bottom of the page you should 
see lines that say:

       Printed in the United States of America
       First Printing

However, it may have the words "Second Printing" 
or "Sixth Printing".

Please let me know which "Printing Number" 
you have, and then I'll know which version of 
the errata I should send to you.

I'll keep in touch Bob,
[-Rick-]

Hi Rick,

> >That's a very good question, but to be honest I haven't considered the
> >time domain characteristics of this function.
>
> Yep, R B-J has a nasty habit of asking tough (but good)
> questions.  Because the freq-domain characteristic is
> rectangular, something like:
>
>    1,2,2,2,...2,2,2,1,0,0,0,0,...,0,0,0,0
>
> I'd guess the time-domain impulse response would
> some sort of complex sequence whose magnitude response
> would be very similar to a sinx(x) sequence.

Based on R B-J's question I took an IFFT of the above sequence, and it
looks as you describe. It looks like the real part is an impulse
followed by 0's and the imaginary part is an asymmetrical sinc() type
sequence.

But if I take the FFT of the above complex sequence and then try to
multiple the resulting frequency domain sequence by the input's
frequency domain sequence, it doesn't look like it is performing a
Hilbert Transformation anymore. This is beyond me, but perhaps
something is lost when the frequency domain coefficients are converted
to time domain and back?

The reason I was doing the above was to see if there is some limit to
the length of the impulse response in the time domain, to get a better
idea of whether or not it can be zero-padded for OLA.

> I'm sorry to say, Bob, but I haven't modeled your
> process of OLA in conjunction with the "freq-domain
> analytic sequence generation' scheme.
> I'll do what I can to model (using MATLAB) this
> whole idea, and get back to you.  Perhaps, together,
> we'll figure this all out.

No problem, and I'm not absolutely certain that OLA is needed. It
appears to be necessary though.

> By the way Bob,
> I'll be happy to send you an errata for
> the 2nd edition of my book if you can tell
> me exactly what is the "Printing Number" of your
> copy of the book.  You can find your
> "Printing Number" on the page just before the
> "Dedication" page.

Thanks, Rick. But there may be no need to send it to me as I bothered
you via email about another issue (Chebyshev window equation) about a
year and a half ago and you sent me the errata then. :) But in case
things have been updated since then, mine is the 3rd printing.

Thanks to both of you for your time and your help.

Regards,
Bob

BobM wrote:
> >
> >    1,2,2,2,...2,2,2,1,0,0,0,0,...,0,0,0,0
> >
> > I'd guess the time-domain impulse response would
> > some sort of complex sequence whose magnitude response
> > would be very similar to a sinx(x) sequence.
>
> Based on R B-J's question I took an IFFT of the above sequence, and it
> looks as you describe.

so what you need to do is window it or something so that it is of
length shorter than your FFT.  that's what you *must* do to use the FFT
and OLA to do convolution.

...
>
> I'm not absolutely certain that OLA is needed. It
> appears to be necessary though.

you can do direct convolution to implement a (delayed) Hilbert
Transformer.


but if you want to do "fast convolution" (by use of the FFT), you need
to insure that your circular convolution does not straddle this nasty
boundary between x[N-1] and x[0] where both are non-zero.  this is why
zero padding is needed. at least for OLA, if it's overlap-save (OLS),
there is another way to do it so that you only save the samples that
weren't contaminated by the bad wrap-around or straddling.  either way
(OLS or OLA) you have to have your frequency response of the FIR filter
be such that the impulse response attached to it is shorter than the
FFT length by enough that you can get some decent quantity of output
samples for each overlapping FFT frame.

> Thanks to both of you for your time and your help.

FWIW.

On 16 Jan 2007 14:00:26 -0800, "BobM" <BobM.DSP@gmail.com> wrote:

>Hi Rick,
>
>> >That's a very good question, but to be honest I haven't considered the
>> >time domain characteristics of this function.
>>
>> Yep, R B-J has a nasty habit of asking tough (but good)
>> questions.  Because the freq-domain characteristic is
>> rectangular, something like:
>>
>>    1,2,2,2,...2,2,2,1,0,0,0,0,...,0,0,0,0
>>
>> I'd guess the time-domain impulse response would
>> some sort of complex sequence whose magnitude response
>> would be very similar to a sinx(x) sequence.
>
>Based on R B-J's question I took an IFFT of the above sequence, and it
>looks as you describe. It looks like the real part is an impulse
>followed by 0's and the imaginary part is an asymmetrical sinc() type
>sequence.
>
>But if I take the FFT of the above complex sequence and then try to
>multiple the resulting frequency domain sequence by the input's
>frequency domain sequence, it doesn't look like it is performing a
>Hilbert Transformation anymore. This is beyond me, but perhaps
>something is lost when the frequency domain coefficients are converted
>to time domain and back?
>
>The reason I was doing the above was to see if there is some limit to
>the length of the impulse response in the time domain, to get a better
>idea of whether or not it can be zero-padded for OLA.
>
>> I'm sorry to say, Bob, but I haven't modeled your
>> process of OLA in conjunction with the "freq-domain
>> analytic sequence generation' scheme.
>> I'll do what I can to model (using MATLAB) this
>> whole idea, and get back to you.  Perhaps, together,
>> we'll figure this all out.

Hi Bob,

>No problem, and I'm not absolutely certain that OLA is needed. It
>appears to be necessary though.

Well, if your time domain sequence is super-long 
then I believe overlap-and-add (OA) or 
overlap-and-save (OS) procesing will be the right 
thing to do.

(We know the freq-domain analytic signal generation 
scheme works, it's just the process of breaking your 
long time sequence into blocks and processing 
each block with a small FFT to implement the "freq-domain 
analytic sig generation" scheme that I'm havin' trouble 
with.)

I have some MATLAB code that implements OS that I wrote.
(That's because I've been convinced by our DSP pal 
Mark Borgerding that OS is more efficient than OA.) 
Anyway, using my OS software and trying implement 
the "freq-domain analytic sig generation" scheme 
on contiguous blocks of time samples has yielded 
results that sound similar to your results.

The real part of my desired analytic signal seems 
like it might be correct, but the imaginary 
part of my time sequence exhibits gross 
discontinuities.  My gut tells me that the OS 
process should work, but I've done a bit of experimenting 
and darn-it I can't get the process to work.
This "bugs" me to no end.  I've been scratching 
my head for some time today, but that hasn't helped 
any.
 
Bob, if we can't figure this out, we can always 
"wave the white flag" and post another message here 
asking for help.  Surely you and I aren't the first 
guys who've tried to solve this problem.  It 
seems so "straightforward"!!

I mean, come on ..., people use FIR filters to generate 
analytic sequences from real sequences.  And OS is an 
efficient way to implement FIR filters.  So how do we 
use OS to generate analytic sequences???
(Shame on me for not knowing this answer.)

There's some subtle "2nd-order effect" taking place here 
that I just can't put my finger on.  Some subtle 
condition, that must be satisfied in our processing,
that we're *not* satisfying. 

>> By the way Bob,
>> I'll be happy to send you an errata for
>> the 2nd edition of my book if you can tell
>> me exactly what is the "Printing Number" of your
>> copy of the book.  You can find your
>> "Printing Number" on the page just before the
>> "Dedication" page.
>
>Thanks, Rick. But there may be no need to send it to me as I bothered
>you via email about another issue (Chebyshev window equation) about a
>year and a half ago and you sent me the errata then. :) But in case
>things have been updated since then, mine is the 3rd printing.

The errata has been updated since we last exchanged 
E-mails.
Bob, do you still have your "...@yahoo.com" E-mail 
address that you used back in May of 2005?

>Thanks to both of you for your time and your help.

Well, I sure wish I had a quick answer for you, 
and I'm sorry that I don't.  I'll keep plugging 
away at this, as time permits.

See Ya',
[-Rick-]

BobM wrote:
> I've been trying to generate an analytic signal (Hilbert
> transformation) in the frequency domain using the method outlined in
> the IEEE Marple paper in the Transactions on Signal Processing
> "Computing the Discrete-Time Analytic Signal via FFT" (9/1999).
...
> I notice that if I use this method for a long input (say several
> seconds of audio) and take an FFT of the entire sample, it seems to
> work except for oddities at the ends.

This method is similar to a rectangular filter in the frequency
domain, which represents a sinc impulse in the time domain.
A sinc is infinite in extent, but the lobes (width and height
envelope) scale inversely with the rectangle width.  So your
long input was a wide rectangle, giving a narrow sinc with
fast decay, leaving only "oddities at the ends".  A short block
input results in a narrower rectangle producing a wider sinc
without as much decay, resulting in the filter impulse wrapping
around your fft aperture and munging the opposite end, perhaps
enough of the ends that there is no unmunged center portion
remaining.

> The real part of the resulting y is always correct. However, the
> imaginary part is only correct if the real input x is a sinusoid with
> either an integer period within the block or an integer period plus
> half of a cycle. For example, the following sinusoids work correctly
> when a block size of 64 is used and i is the index:
>
> x=sin(2*pi*(1/64)*i)
> x=sin(2*pi*(1.5/64)*i)

Note that these functions are continuous with the zero padding
(they end with x=0 and so does the zero padding) in your window.


IMHO. YMMV.
-- 
rhn A.T nicholson d.0.t C-o-M

Rick Lyons wrote:
> (That's because I've been convinced by our DSP pal
> Mark Borgerding that OS is more efficient than OA.)

i dunno why it wouldn't.  OS and OA have the same computation for the
FFT, applying frequency response, and iFFT.  and they both advance your
input pointer the same amount per frame.  and the OA requires padding
the input frame with zeros and adding the overlapped output frames.  OS
requires copying some samples over from the previous frame instead of
padding and then just picking the segment of samples that is not
contaminated from the impulse response straddling the wrap around point
between x[N-1] and x[0].  copying samples is maybe a little slower than
stuffing zero but i would not expect the cost of those extra reads to
be as much as a the same number of reads, adds, and write.  actually OA
would be twice the number of reads (twice the overlap length) than OS.
 i'm convinced, too (but i always thunk it).  i wonder if, for
numerical reasons with finite precision, if OA might "sound better".
the "discontinuity" between x[N-1] and x[0] gets in the math in OS and
sorta comes out in the wash (except for those L-1 samples that are not
saved), but i wonder if there is some numerical "lingering" of that
"click" if the word size is small enough even in the samples you keep.

> Anyway, using my OS software and trying implement
> the "freq-domain analytic sig generation" scheme
> on contiguous blocks of time samples has yielded
> results that sound similar to your results.

even if the impulse response is no longer than N-M+1 where N is the FFT
length and M is your hop size (the number of samples you keep each
frame)? 

r b-j