john wrote: ...> One way to split this real signal into quadrature channels is to apply > two bandpass filters constructed from lowpass prototypes. Make a FIR > LPF with half your desired bandwidth. Then multiply the taps by sine > and cosine of 240 kHz to make them into BPFs. If you pass the input > through each of the filters (sine and cosine), the outputs will be in > quadrature.I think you left out a crucial step. One takes sines and cosines of angles, not frequencies. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Hilbert transform using FFT approach
Started by ●September 21, 2005
Reply by ●September 23, 20052005-09-23
Reply by ●September 23, 20052005-09-23
Jerry Avins wrote:> john wrote: > > ... > > > One way to split this real signal into quadrature channels is to apply > > two bandpass filters constructed from lowpass prototypes. Make a FIR > > LPF with half your desired bandwidth. Then multiply the taps by sine > > and cosine of 240 kHz to make them into BPFs. If you pass the input > > through each of the filters (sine and cosine), the outputs will be in > > quadrature. > > I think you left out a crucial step. One takes sines and cosines of > angles, not frequencies. > > Jerry > -- > Engineering is the art of making what you want from things you can get. > =AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF==AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF= =AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF=AF Picky picky. Here is some Matlab code: Fc=3D240e3;Fs=3D10e6;W=3D100e3; x=3Dsin(2*pi*[0:999]*Fc/Fs); Ntaps=3D200; hlp=3Dfir1(Ntaps-1,(W/2)/(Fs/2)); hbps=3D2*hlp.*sin(2*pi*[0:Ntaps-1]*Fc/Fs); hbpc=3D2*hlp.*cos(2*pi*[0:Ntaps-1]*Fc/Fs); x =3D sin(2*pi*(Fc-W/4)*[0:999]/Fs); xs =3D filter(hbps,1,x); xc =3D filter(hbpc,1,x); plot([xs' xc']);grid;shg; John
Reply by ●September 23, 20052005-09-23
john wrote: ...> Picky picky. Here is some Matlab code:I wasn't trying to be picky. Without that code, I would have had to guess or look it up. (I would have guessed right. I have the means to look.) You were trying to teach an even less informed guy than I. Filling in the details was important. ... Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●September 27, 20052005-09-27
On 23 Sep 2005 10:02:49 -0700, "john" <johns@xetron.com> wrote:> >Jerry Avins wrote: >> john wrote: >> >> ... >> >> > One way to split this real signal into quadrature channels is to apply >> > two bandpass filters constructed from lowpass prototypes. Make a FIR >> > LPF with half your desired bandwidth. Then multiply the taps by sine >> > and cosine of 240 kHz to make them into BPFs. If you pass the input >> > through each of the filters (sine and cosine), the outputs will be in >> > quadrature. >> >> I think you left out a crucial step. One takes sines and cosines of >> angles, not frequencies. >> >> Jerry> >Picky picky. Here is some Matlab code: > >Fc=3D240e3;Fs=3D10e6;W=3D100e3; >x=3Dsin(2*pi*[0:999]*Fc/Fs); >Ntaps=3D200; >hlp=3Dfir1(Ntaps-1,(W/2)/(Fs/2)); >hbps=3D2*hlp.*sin(2*pi*[0:Ntaps-1]*Fc/Fs); >hbpc=3D2*hlp.*cos(2*pi*[0:Ntaps-1]*Fc/Fs); >x =3D sin(2*pi*(Fc-W/4)*[0:999]/Fs); >xs =3D filter(hbps,1,x); >xc =3D filter(hbpc,1,x); >plot([xs' xc']);grid;shg; > >JohnHi John, I don't think Jerry was being picky. He was just being correct. Ya' know what I found out about that "multiply a lowpass filter by sines & cosines" method for designing a Hilbert transformer? That method sometimes gives you bandpass filters whose phase responses are not exactly linear. After removing those strange "3D" characters in your code, I ran your code. Your example is a good example of the point you were making. It might be interesting for you to experiment with far fewer than 200 taps and different BPF center frequencies to see the phase nonlinearity. For example: Fc=1.5e6; Fs=10e6; W=100e3; x=sin(2*pi*[0:999]*Fc/Fs); Ntaps=17; hlp=fir1(Ntaps-1,(W/2)/(Fs/2)); hbps=2*hlp.*sin(2*pi*[0:Ntaps-1]*Fc/Fs); hbpc=2*hlp.*cos(2*pi*[0:Ntaps-1]*Fc/Fs); x = sin(2*pi*(Fc-W/4)*[0:999]/Fs); xs = filter(hbps,1,x); xc = filter(hbpc,1,x); figure(1) freqz(hbps,1024); grid on, zoom on figure(2) subplot(2,1,1), plot(hbps,'-bo'), grid on, zoom on subplot(2,1,2), plot(grpdelay(hbps,1,1024)); grid on, zoom on figure(3) plot([xs' xc']); grid on, zoom on I guess using hundreds of taps, as you did, drastically reduces phase nonlinearity problems. Interesting, ... huh? Thanks for making me think about this again. See Ya', [-Rick-]
Reply by ●September 27, 20052005-09-27
On Wed, 21 Sep 2005 18:36:56 GMT, Al Clark <dsp@danvillesignal.com> wrote: (snipped)>You can use an FFT to generate Hilbert transforms but it is often not the >easiest way. > >The way we usually generate a Hilbert transform is two use >antisymmetrical odd length FIR filters. This method is clearly discussed >in Rick Lyon's book: Understanding Digital Signal Processing (second >edition). > >We implemented a hilbert tranformer for essentially this same purpose >last week for an AM demodulator application running on one of our SHARC >based DSP boards. You often need a square root for envelope detectors. >This is very easy with a SHARC since it has an instruction for this >purpose. >Hi Al, you sure piqued my interest in learning about this "Sharc instruction" that's useful for envelope detection. I went to the Analog Dev. website to see what I could learn, but there's SOOoo much stuff on their website, that I was unsuccessful. (By "unsuccessful" I mean that I got sidetracked reading their "white papers" and "technical articles". Ha ha.) So Al, I'd like to learn the details of that "Sharc instruction" and see exactly how it works (if that's possible). Is there a website where I can learn about that instruction? Thanks Al, [-Rick-]
Reply by ●September 27, 20052005-09-27
Rick Lyons wrote:> On 23 Sep 2005 10:02:49 -0700, "john" <johns@xetron.com> wrote: > > > > >Jerry Avins wrote: > >> john wrote: > >> > >> ... > >> > >> > One way to split this real signal into quadrature channels is to apply > >> > two bandpass filters constructed from lowpass prototypes. Make a FIR > >> > LPF with half your desired bandwidth. Then multiply the taps by sine > >> > and cosine of 240 kHz to make them into BPFs. If you pass the input > >> > through each of the filters (sine and cosine), the outputs will be in > >> > quadrature. > >> > >> I think you left out a crucial step. One takes sines and cosines of > >> angles, not frequencies. > >> > >> Jerry > > > > >Picky picky. Here is some Matlab code: > > > >Fc=3D240e3;Fs=3D10e6;W=3D100e3; > >x=3Dsin(2*pi*[0:999]*Fc/Fs); > >Ntaps=3D200; > >hlp=3Dfir1(Ntaps-1,(W/2)/(Fs/2)); > >hbps=3D2*hlp.*sin(2*pi*[0:Ntaps-1]*Fc/Fs); > >hbpc=3D2*hlp.*cos(2*pi*[0:Ntaps-1]*Fc/Fs); > >x =3D sin(2*pi*(Fc-W/4)*[0:999]/Fs); > >xs =3D filter(hbps,1,x); > >xc =3D filter(hbpc,1,x); > >plot([xs' xc']);grid;shg; > > > >John > > Hi John, > > I don't think Jerry was being picky. > He was just being correct. > > Ya' know what I found out about that > "multiply a lowpass filter by sines & cosines" > method for designing a Hilbert transformer? > That method sometimes gives you bandpass > filters whose phase responses > are not exactly linear. > > After removing those strange "3D" characters in your > code, I ran your code. Your example is a good > example of the point you were making. > > It might be interesting > for you to experiment with far fewer than 200 taps > and different BPF center frequencies to > see the phase nonlinearity. For example: > > Fc=1.5e6; > Fs=10e6; > W=100e3; > x=sin(2*pi*[0:999]*Fc/Fs); > Ntaps=17; > hlp=fir1(Ntaps-1,(W/2)/(Fs/2)); > hbps=2*hlp.*sin(2*pi*[0:Ntaps-1]*Fc/Fs); > hbpc=2*hlp.*cos(2*pi*[0:Ntaps-1]*Fc/Fs); > x = sin(2*pi*(Fc-W/4)*[0:999]/Fs); > xs = filter(hbps,1,x); > xc = filter(hbpc,1,x); > > figure(1) > freqz(hbps,1024); grid on, zoom on > figure(2) > subplot(2,1,1), plot(hbps,'-bo'), grid on, zoom on > subplot(2,1,2), plot(grpdelay(hbps,1,1024)); grid on, zoom on > figure(3) > plot([xs' xc']); grid on, zoom on > > I guess using hundreds of taps, as you did, > drastically reduces phase nonlinearity problems. > Interesting, ... huh? > > Thanks for making me think about this again. > > See Ya', > [-Rick-]Yes, I see. I wonder if the nonlinear phase occurs when Fs/Fc is not an integer? John
Reply by ●September 27, 20052005-09-27
R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote in news:433914ac.523428515@news.sf.sbcglobal.net:> On Wed, 21 Sep 2005 18:36:56 GMT, Al Clark <dsp@danvillesignal.com> > wrote: > > (snipped) > >>You can use an FFT to generate Hilbert transforms but it is often not >>the easiest way. >> >>The way we usually generate a Hilbert transform is two use >>antisymmetrical odd length FIR filters. This method is clearly >>discussed in Rick Lyon's book: Understanding Digital Signal Processing >>(second edition). >> >>We implemented a hilbert tranformer for essentially this same purpose >>last week for an AM demodulator application running on one of our >>SHARC based DSP boards. You often need a square root for envelope >>detectors. This is very easy with a SHARC since it has an instruction >>for this purpose. >> > > Hi Al, > > you sure piqued my interest in learning about this > "Sharc instruction" that's useful for envelope detection. > > I went to the Analog Dev. website to see what I could > learn, but there's SOOoo much stuff on their > website, that I was unsuccessful. > > (By "unsuccessful" I mean that I got sidetracked > reading their "white papers" and "technical > articles". Ha ha.) > > So Al, I'd like to learn the details of that "Sharc instruction" > and see exactly how it works (if that's possible). > > Is there a website where I can learn about that instruction? > > Thanks Al, > [-Rick-] > >The envelope is SQRT (I^2 + R^2) 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. -- Al Clark Danville Signal Processing, Inc. -------------------------------------------------------------------- Purveyors of Fine DSP Hardware and other Cool Stuff Available at http://www.danvillesignal.com
Reply by ●September 27, 20052005-09-27
Al Clark wrote: ...> The envelope is SQRT (I^2 + R^2) > > > 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.IIRC, each iteration of the algorithm doubles the number of good bits. So if RSQRTS gives 4 bits, the first iteration gives 8, the second gives 16, and the third gives 32. If the result of RSQRTS is only a little better, 32 could become 40. I assume the number of iterations is controlled by a loop counter, so that the instruction count doesn't depend on how quickly the calculation converges. Is that right? Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●September 27, 20052005-09-27
"Jerry Avins" <jya@ieee.org> wrote in message news:WO2dnR0l_OIZx6TeRVn-oA@rcn.net...> Al Clark wrote: > > ... > >> The envelope is SQRT (I^2 + R^2) >> >> >> 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. > > IIRC, each iteration of the algorithm doubles the number of good bits. So if > RSQRTS gives 4 bits, the first iteration gives 8, the second gives 16, and the > third gives 32. If the result of RSQRTS is only a little better, 32 could > become 40. > > I assume the number of iterations is controlled by a loop counter, so that the > instruction count doesn't depend on how quickly the calculation converges. Is > that right?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.
Reply by ●September 28, 20052005-09-28
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? > >JohnHi, 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-]
