Hilbert transform using FFT approach

On Tue, 27 Sep 2005 14:01:33 GMT, Al Clark <dsp@danvillesignal.com>
wrote:

  (snipped)

Hi Al, 

>The envelope is SQRT (I^2 + R^2)

Yep.

>The instruction is fn = RSQRTS fx
>
>The instruction creates a 4 bit accurate floating point seed for 
>SQRT(1/fx). This is the starting point for a calculation using the 
>Newton-Raphson iteration algorithm. You get the SQRT by multiplying fn * 
>fx.
>
>The SHARC Instruction manual has example code for this algorithm.
>
>Since the seed is of known accuracy, you know the result will be accurate  
>after a specific number of iterations. It takes 13 instructions to get to 
>32 or 40 bit accuracy. 

Ah ha.  OK.  Yes. Interesting, very interesting.

Thanks very much Al,
[-Rick-]

On Tue, 27 Sep 2005 15:37:13 GMT, "Jon Harris"
<jon99_harris7@hotmail.com> wrote:

  (snipped)
>
>In the SHARC example code, there is no loop counter.  Instead, it is all 
>"unrolled" code.  Doing this saves the 1 instruction required to set-up the 
>loop, which could fairly significant.  There is a note that says "leave out 
>these instructions" for a reduced accuracy version.
>
>Rick, see 
>http://www.analog.com/UploadedFiles/Associated_Docs/3307288358013isr_21xxx_07.pdf, 
>p.46-47 for complete details. 

Hi Jon,

   Hey, thanks much for the URL!

I see the three iterations of the Newton-Raphson 
Inverse method,

p(n+1) = 0.5*p(n)[3 - x*p(n)^2],

in their code comments.    Neat!

As far as I can tell, they use use two iterations 
for 32-bit floating-pt data and three iterations for 
40-bit floating-pt data.  In both cases their error is 
no greater than +- 1 LSB.

Thanks again for that URL Jon.
That PDF file was interesting to scan through.
(I haven't looked at an assembly intruction 
for years and years.)

See Ya',
[-Rick-]

R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote in
news:433a5863.606298625@news.sf.sbcglobal.net: 

> On Tue, 27 Sep 2005 15:37:13 GMT, "Jon Harris"
><jon99_harris7@hotmail.com> wrote:
> 
>   (snipped)
>>
>>In the SHARC example code, there is no loop counter.  Instead, it is
>>all "unrolled" code.  Doing this saves the 1 instruction required to
>>set-up the loop, which could fairly significant.  There is a note that
>>says "leave out these instructions" for a reduced accuracy version.
>>
>>Rick, see 
>>http://www.analog.com/UploadedFiles/Associated_Docs/3307288358013isr_21
>>xxx_07.pdf, p.46-47 for complete details. 
> 
> Hi Jon,
> 
>    Hey, thanks much for the URL!
> 
> I see the three iterations of the Newton-Raphson 
> Inverse method,
> 
> p(n+1) = 0.5*p(n)[3 - x*p(n)^2],
> 
> in their code comments.    Neat!
> 
> As far as I can tell, they use use two iterations 
> for 32-bit floating-pt data and three iterations for 
> 40-bit floating-pt data.  In both cases their error is 
> no greater than +- 1 LSB.
> 
> Thanks again for that URL Jon.
> That PDF file was interesting to scan through.
> (I haven't looked at an assembly intruction 
> for years and years.)

Time to start? Who says you can't teach an old dsp dog a new trick?

> 
> See Ya',
> [-Rick-]
> 
> 



-- 
Al Clark
Danville Signal Processing, Inc.
--------------------------------------------------------------------
Purveyors of Fine DSP Hardware and other Cool Stuff
Available at http://www.danvillesignal.com

"Rick Lyons" <R.Lyons@_BOGUS_ieee.org> wrote in message 
news:433a5863.606298625@news.sf.sbcglobal.net...
> On Tue, 27 Sep 2005 15:37:13 GMT, "Jon Harris"
> <jon99_harris7@hotmail.com> wrote:
>
>  (snipped)
>>
>>In the SHARC example code, there is no loop counter.  Instead, it is all
>>"unrolled" code.  Doing this saves the 1 instruction required to set-up the
>>loop, which could fairly significant.  There is a note that says "leave out
>>these instructions" for a reduced accuracy version.
>>
>>Rick, see
>>http://www.analog.com/UploadedFiles/Associated_Docs/3307288358013isr_21xxx_07.pdf,
>>p.46-47 for complete details.
>
> Hi Jon,
>
>   Hey, thanks much for the URL!
>
> I see the three iterations of the Newton-Raphson
> Inverse method,
>
> p(n+1) = 0.5*p(n)[3 - x*p(n)^2],
>
> in their code comments.    Neat!
>
> As far as I can tell, they use use two iterations
> for 32-bit floating-pt data and three iterations for
> 40-bit floating-pt data.  In both cases their error is
> no greater than +- 1 LSB.

Yes, that is exactly right.  The whole key to the algorithm is coming up with a 
"good" starting point, which the hardware-based look-up table does nicely (and I 
assume with a fairly small look-up table).  FYI, they have a similar instruction 
and algorithm for computing 1/X (RECIPS), which allows for fairly quick division 
routines.  It is described on p. 44-45 of the same document.  Those two little 
instructions are some of my favorite things about the SHARC!

> Thanks again for that URL Jon.
> That PDF file was interesting to scan through.
> (I haven't looked at an assembly intruction
> for years and years.)

You're welcome. 

Jon Harris wrote:

   ...

>        The whole key to the algorithm is coming up with a 
> "good" starting point, which the hardware-based look-up table does nicely (and I 
> assume with a fairly small look-up table).  FYI, they have a similar instruction 
> and algorithm for computing 1/X (RECIPS), which allows for fairly quick division 
> routines.  It is described on p. 44-45 of the same document.  Those two little 
> instructions are some of my favorite things about the SHARC!

With floating-point numbers, halving the exponent is a good start on the 
seed. With non-offset formats, that's just a shift. The carry out of 
that shift and the two MSBs can index into an 8-element LUT that gives a 
5-bit sqrt mantissa. I guess from that that SHARC's LUT has 4 entries.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

Hi Rick,

This may have already been said. I am getting in late and am not going
to read all 30 posts. If so, I appologize.

If you start with a symmetric LP filter and delay the sine and cosine
by the delay of the lowpass filter before using them to multiply the LP
filter then the resulting filters will be symmetric/antisymmetric and
have linear phase guaranteed. Properly aligned windowing at any point
in the process won't change that.

Does that take care of the problem?

Dirk

Rick,

Looking at the signal (center?) freq and the sample rate, it seems to
be begging for a really efficient multirate method (mix, decimate,
hilbert transform, interpolate, remix/make real).

For envelope detection only, maybe you could even do the sqrt at the
decimated rate, do a real interpolation, and skip the remix/make real
step.  Even more efficient.

Dirk

Rick Lyons wrote:
> On 27 Sep 2005 03:59:35 -0700, "john" <johns@xetron.com> wrote:
>
> >
>   (snipped)
> >
> >Yes, I see. I wonder if the nonlinear phase occurs when Fs/Fc is not an
> >integer?
> >
> >John
>
> Hi,
>
> Oh shoot.  Such a good question!!
>
> I don't have a quick answer off the top of my
> head.  I just remember that multiplying a sequence
> of symmetrical LP filter coefficients by a cos & sin
> sequence could yield bandpass coeffs that were not
> symmetrical.
>
> Instead of using the words "could yield" in the above
> sentence, I propbably should have used the phrase
> "almost always yields".
>
> I'll bet a bottle of beer that (at least) two things
> affect the BP coeffs' linearity:
>
> 1) Fs/Fc ratio (just as you say)
> 2) the length of the original lowpass coeff sequence.
>
> Ah, I just now dug up an old note that I had written.
> It says that if Fs/Fc = 4, and the number of
> lowpass filter taps is odd, then the bandpass
> coeffs will be symmetrical (linear phase).
>
> Of course, windowing nonlinear BP filter coeffs will
> improve their linearity.
>
> Hey John, I just had a (crazy) thought.  I wonder if
> we could design a slightly nonlinear phase lowpass
> filter and then have the frequency translation operation
> cancel that nonlinearity to provide symmetrical
> bandpass filter coefficients?
>
> This sounds like a challenging (however impractical)
> homework problem.
> 
> See Ya',
> [-Rick-]

On 28 Sep 2005 08:50:02 -0700, "dbell" <dbell@niitek.com> wrote:

>Hi Rick,

Hi Dirk,

>This may have already been said. I am getting in late and am not going
>to read all 30 posts. If so, I appologize.

   Ha.  I haven't read all the post either.
(I'm just thinking about turning a lowpass filter
into a bandpass filter.)

>If you start with a symmetric LP filter and delay the sine and cosine
>by the delay of the lowpass filter before using them to multiply the LP
>filter then the resulting filters will be symmetric/antisymmetric and
>have linear phase guaranteed. Properly aligned windowing at any point
>in the process won't change that.

Humm, as that always true?  Here's my thinking: if 
K samples of a sinewave sequence are not symmetrical 
then K samples of that sinewave delayed in time will 
also not be symmetrical.

Am I missing something here? 

Yes.  The alignment between the sine and cosine must of course be
fixed, but how you position them against the lowpass filter taps before
you multiply is up to you.  You can make the resulting filters
symmetrical by effectively shifting the sine positive zero crossing and
cosine maximum to both align with the point of symmetry of the original
symmetric lowpass filter. 

Dirk

On 28 Sep 2005 14:19:09 -0700, "dbell" <dbell@niitek.com> wrote:

>Yes.  The alignment between the sine and cosine must of course be
>fixed, but how you position them against the lowpass filter taps before
>you multiply is up to you.  You can make the resulting filters
>symmetrical by effectively shifting the sine positive zero crossing and
>cosine maximum to both align with the point of symmetry of the original
>symmetric lowpass filter. 
>
>Dirk

Hi Dirk,

  Humm, maybe am missing something, and maybe I'm not.

Let's just talk about a cosine sequence.  I was 
thinking about a cosine sequence of some 
frequency Fc sampled at a sample rate of Fs.  What 
I was thinking was that this sequence will not 
be symmetrical for all possible sequence lengths N.

Dirk, I'll be out of town for the next 7-10 days on business 
so if you reply and don't hear back from me for a while, 
please know that I'm not ignoring you.

[-Rick-]