DSPRelated.com
Forums

Generating digital chirp signal

Started by Unknown September 23, 2007
Hi!

I'm trying to implement a digital chirp signal in hardware.
The formula is:

c[n] = exp(j*ang[n]), n=0,1,2,...,N-1
ang[n] = (pi*n^2) / N.

N is any *prime* number such that N < 1320.

I have an idea to do that: sampling the exp(j*k) ant insert the
samples into a look up table,
I'll need to calculate an ang[n] every time so c[n] = LUT{ang[n]}.
(LUT = Look UP Table).

But I don't like this idea because it's very generic and not
efficient.

I thought about using "chirp z-transform" to generate the digital
signal chirp, but
I don't how...

Any suggestions? Any other ideas to generate the signal?

Thanks in advance

Hi Benome,
What you are proposing can be quite efficient and is similar to the 
standard technique used for Direct Digital Synthesis, except that you 
need a progressibely increasing phase-increment instead of a constant 
phase-increment.

To generate the phase-increment, you could generate successive n^2 
values and then multiply by each by pi/N as you suggest, but in fact 
there is a simpler way which uses second differences and requires no 
multiplies.   Actually, you don't need to generate n^2 explicitly at 
all, as you will see in the following routine:-

CHIRP ROUTINE
Initially:   ang =0,   D1 = -pi/n,   D2 = 2*pi/N
These variables are all fixed-point, with an integer part and a 
fractional part.
   1.  Look up cos(ang) and sin(ang)to give you the required exp(j*ang).
   2.  D1 = D1 + D2
   3.  ang = ang + D1
   4.  If ang < ang(max), go to (1)
   5.  End of chirp.

The angle (ang) needs to be held in a long register.  The high-order 
bits are used to access the LUT and the low-order bits are there to 
minimise the accumulated phase error.


By the way, you can reduce the size of the LUT by a factor of eight by 
taking advantage of the symmetry of the cos and sin functions.  Have a 
LUT for only cos(ang), and then only in the range to 0 to pi/2.  You 
then use appropriate rules to find cos(arg) and sin(arg) for any value 
of (arg).

Regards,
John



benomed@gmail.com wrote:
> Hi! > > I'm trying to implement a digital chirp signal in hardware. > The formula is: > > c[n] = exp(j*ang[n]), n=0,1,2,...,N-1 > ang[n] = (pi*n^2) / N. > > N is any *prime* number such that N < 1320. > > I have an idea to do that: sampling the exp(j*k) ant insert the > samples into a look up table, > I'll need to calculate an ang[n] every time so c[n] = LUT{ang[n]}. > (LUT = Look UP Table). > > But I don't like this idea because it's very generic and not > efficient. > > I thought about using "chirp z-transform" to generate the digital > signal chirp, but > I don't how... > > Any suggestions? Any other ideas to generate the signal? > > Thanks in advance >
On Sep 23, 5:58 pm, John Monro <johnmo...@optusnet.removethis.com.au>
wrote:
> Hi Benome, > What you are proposing can be quite efficient and is similar to the > standard technique used for Direct Digital Synthesis, except that you > need a progressibely increasing phase-increment instead of a constant > phase-increment.
but that's the "LUT" method and he said he didn't like it. (was interpolation the reason the OP didn't like it?)
> To generate the phase-increment, you could generate successive n^2 > values and then multiply by each by pi/N as you suggest, but in fact > there is a simpler way which uses second differences and requires no > multiplies. Actually, you don't need to generate n^2 explicitly at > all, as you will see in the following routine:- > > CHIRP ROUTINE > Initially: ang =0, D1 = -pi/n, D2 = 2*pi/N > These variables are all fixed-point, with an integer part and a > fractional part. > 1. Look up cos(ang) and sin(ang)to give you the required exp(j*ang). > 2. D1 = D1 + D2 > 3. ang = ang + D1 > 4. If ang < ang(max), go to (1) > 5. End of chirp. > > The angle (ang) needs to be held in a long register. The high-order > bits are used to access the LUT and the low-order bits are there to > minimise the accumulated phase error. > > By the way, you can reduce the size of the LUT by a factor of eight by > taking advantage of the symmetry of the cos and sin functions. Have a > LUT for only cos(ang), and then only in the range to 0 to pi/2. You > then use appropriate rules to find cos(arg) and sin(arg) for any value > of (arg).
one thing to think about is that the chirp signal is not bandlimited (even if the range of the instantaneous frequency of the chirp always had magnitude less than Nyquist). sampling it (which is what you're doing if your interpolation gave you the precise time-domain values of the sin(t^2) function. theoretically just outputting those LUT values, even if the bandlimited interpolation was perfect (using an arbitrarily long sinc(.) function), it's not the same as sampling a "properly" bandlimited chirp. you might have to think about images and aliases if the chirp sweep rate is fast enough. r b-j
robert bristow-johnson wrote:
> On Sep 23, 5:58 pm, John Monro <johnmo...@optusnet.removethis.com.au> > wrote: >> Hi Benome, >> What you are proposing can be quite efficient and is similar to the >> standard technique used for Direct Digital Synthesis, except that you >> need a progressibely increasing phase-increment instead of a constant >> phase-increment. > > but that's the "LUT" method and he said he didn't like it. (was > interpolation the reason the OP didn't like it?)
Dunno! The Op was looking for suggestions. The fact that I reduced the size of the LUT by a factor of eight may make him reconsider using a LUT. I don't quite understand your references to interpolation, as the method I suggested does not use interpolation.
> >> To generate the phase-increment, you could generate successive n^2 >> values and then multiply by each by pi/N as you suggest, but in fact >> there is a simpler way which uses second differences and requires no >> multiplies. Actually, you don't need to generate n^2 explicitly at >> all, as you will see in the following routine:- >> >> CHIRP ROUTINE >> Initially: ang =0, D1 = -pi/n, D2 = 2*pi/N >> These variables are all fixed-point, with an integer part and a >> fractional part. >> 1. Look up cos(ang) and sin(ang)to give you the required exp(j*ang). >> 2. D1 = D1 + D2 >> 3. ang = ang + D1 >> 4. If ang < ang(max), go to (1) >> 5. End of chirp. >> >> The angle (ang) needs to be held in a long register. The high-order >> bits are used to access the LUT and the low-order bits are there to >> minimise the accumulated phase error. >> >> By the way, you can reduce the size of the LUT by a factor of eight by >> taking advantage of the symmetry of the cos and sin functions. Have a >> LUT for only cos(ang), and then only in the range to 0 to pi/2. You >> then use appropriate rules to find cos(arg) and sin(arg) for any value >> of (arg). > > one thing to think about is that the chirp signal is not bandlimited > (even if the range of the instantaneous frequency of the chirp always > had magnitude less than Nyquist). sampling it (which is what you're > doing if your interpolation gave you the precise time-domain values of > the sin(t^2) function. theoretically just outputting those LUT > values, even if the bandlimited interpolation was perfect (using an > arbitrarily long sinc(.) function), it's not the same as sampling a > "properly" bandlimited chirp. you might have to think about images > and aliases if the chirp sweep rate is fast enough. > > r b-j >
I agree that this chirp signal is not bandlimited, but do not understand the reference to "...your interpolation.." Regards, John
On Sep 24, 4:46 am, John Monro <johnmo...@optusnet.removethis.com.au>
wrote:
> robert bristow-johnson wrote: > > On Sep 23, 5:58 pm, John Monro <johnmo...@optusnet.removethis.com.au> > > wrote: > >> Hi Benome, > >> What you are proposing can be quite efficient and is similar to the > >> standard technique used for Direct Digital Synthesis, except that you > >> need a progressibely increasing phase-increment instead of a constant > >> phase-increment. > > > but that's the "LUT" method and he said he didn't like it. (was > > interpolation the reason the OP didn't like it?) > > Dunno! The Op was looking for suggestions. The fact that I reduced the > size of the LUT by a factor of eight may make him reconsider using a LUT. > > I don't quite understand your references to interpolation, as the method > I suggested does not use interpolation. > > > > > > >> To generate the phase-increment, you could generate successive n^2 > >> values and then multiply by each by pi/N as you suggest, but in fact > >> there is a simpler way which uses second differences and requires no > >> multiplies. Actually, you don't need to generate n^2 explicitly at > >> all, as you will see in the following routine:- > > >> CHIRP ROUTINE > >> Initially: ang =0, D1 = -pi/n, D2 = 2*pi/N > >> These variables are all fixed-point, with an integer part and a > >> fractional part. > >> 1. Look up cos(ang) and sin(ang)to give you the required exp(j*ang). > >> 2. D1 = D1 + D2 > >> 3. ang = ang + D1 > >> 4. If ang < ang(max), go to (1) > >> 5. End of chirp. > > >> The angle (ang) needs to be held in a long register. The high-order > >> bits are used to access the LUT and the low-order bits are there to > >> minimise the accumulated phase error. > > >> By the way, you can reduce the size of the LUT by a factor of eight by > >> taking advantage of the symmetry of the cos and sin functions. Have a > >> LUT for only cos(ang), and then only in the range to 0 to pi/2. You > >> then use appropriate rules to find cos(arg) and sin(arg) for any value > >> of (arg). > > > one thing to think about is that the chirp signal is not bandlimited > > (even if the range of the instantaneous frequency of the chirp always > > had magnitude less than Nyquist). sampling it (which is what you're > > doing if your interpolation gave you the precise time-domain values of > > the sin(t^2) function. theoretically just outputting those LUT > > values, even if the bandlimited interpolation was perfect (using an > > arbitrarily long sinc(.) function), it's not the same as sampling a > > "properly" bandlimited chirp. you might have to think about images > > and aliases if the chirp sweep rate is fast enough. > > I agree that this chirp signal is not bandlimited, but do not understand > the reference to "...your interpolation.."
in step #1:
> 1. Look up cos(ang) and sin(ang)to give you the required exp(j*ang).
how do you "look up" when ang takes on a value that, when scaled to be an index into the LUT, is not precisely an integer value (so a single look-up isn't precise)? what do you do with "the low-order bits [that] are there to minimise the accumulated phase error"? i presume that step #1 has some interpolation intrinsic to the operation to handle such values of "ang", since to reference to this is made in the other steps. now, that's okay John, i'm not criticizing. i'm only pointing out that whatever you're doing with these "low-order bits", even simply ignoring them and throwing them away (regarding table lookup, i know that you keep them in the phase accumulator so that you have some decent precision on the instantaneous frequency), whatever reasonable thing you do with those bits in the table lookup (incl. ignoring them) is a form of interpolation. now let's suppose, for the point of illustration is that what you do with those bits is apply them (and the neighboring values in the table) to an idealized reconstruction formula like: +M/2 x(m+a) = SUM{ x[(m+n)modulo_N] * sinc(a-n) } n=-M/2+1 where x[n] are the discrete samples of the sine in the table and for an integer, n, x(n) = x[n] (with no meaningful loss of generality, i'm normalizing T for simplicity) and M is a very large positive integer and is the number of terms of the summation in the interpoation. x(t) is the continuous-time output. m+a is the "ang" parameter where m is the integer part (m = floor(m+a)) or the "high-order bits" and "a" is the fractional part (a = (m+a) - ang(m+a) and 0 <= a < 1). also x[n] are the samples for a single cycle of a sinusoid function: x[n] = cos(2*pi*n/N); or, conceptually, x[n] = e^(j*2*pi*n/N); (j^2 = -1) where N is the size of the lookup table and (m+n)modulo_N = (m+n) - N*floor((m+n)/N) . suppose you do that for interpolation in the LUT. now that original sinusoid represented in the table is clearly bandlimited to much less than Nyquist, so the ideal sinc-based reconstruction should give you pretty much exactly the sinusoid an *any* phase (besides phases that are 2*pi*n/N) and it shouldn't matter at all what the phase is or how it is increasing (like a t^2 function). what should come out of +M/2 x(m+a) = SUM{ x[(m+n)modulo_N] * sinc(a-n) } n=-M/2+1 should be awful damn close to x(m+a) = e^(j*2*pi*(m+a)/N) and, in fact, can come out exact if M=N and we change the sinc(.) to the Dirichlet thingie: +N/2 x(m+a) = SUM{x[(m+n)mod_N] * sin(pi*(a-n))/(N*sin(pi/N*(a-n)))} n=-N/2+1 (i hope Google Groups didn't wrap that nastily.) i'm still debating with myself if, for even N, it should be sin(pi*t)/ (N*sin(pi/N*t)) or sin(pi*t)/(N*tan(pi/N*t)) for exact reconconstruction (and i wish someone can settle this and persuade me of the answer). now, if *that* idealized interpolation is what you are doing for step #1, then we know we are taking this ideal band- UNlimited chirp and sampling it at equal spaced times. but since it's not bandlimited, what are the consequences of sampling it? if your interpolation is less ideal, it's even a more difficult question to how nasty the chirp output will be. maybe not very nasty if N is large and the sweep rate is sufficiently slow. r b-j
On Mon, 24 Sep 2007 06:08:46 -0000, robert bristow-johnson
<rbj@audioimagination.com> wrote:

>On Sep 23, 5:58 pm, John Monro <johnmo...@optusnet.removethis.com.au> >wrote: >> Hi Benome, >> What you are proposing can be quite efficient and is similar to the >> standard technique used for Direct Digital Synthesis, except that you >> need a progressibely increasing phase-increment instead of a constant >> phase-increment. > >but that's the "LUT" method and he said he didn't like it. (was >interpolation the reason the OP didn't like it?)
It doesn't need to be interpreted as a LUT implementation. The key is just that the phase is accumulated quadratically (for the chirp) instead of linearly (for a DDS). One can just use the accumulated phase as the argument to a sin or cos function, and, voila! a chirp. When I started out in radar signal processing this was routine, and where I worked it was called a DDA (Digital Differential Analyzer) instead of a DDS. Eric Jacobsen Minister of Algorithms Abineau Communications http://www.ericjacobsen.org
>> but since it's not bandlimited, what are the consequences of sampling
it? Aliasing. The parts of the spectrum that don't fit into the bandwidth simply wrap around. As was said, the ideal spectrum is infinitely wide. If I put two ramps end-to-end, it looks almost like FM modulation (triangle instead of sine). Therefore, standard FM radio equations might come in handy to estimate the effective bandwidth (for example Carson's rule, note the disclaimer for discontinuities: http://en.wikipedia.org/wiki/Carson_bandwidth_rule). This is only a crude guess, having no suitable (radar?) textbook around. It might be that also the interpolation error from the wavetable appears in a similar way (as phase modulation) around the instantaneous frequency. -mn
Eric Jacobsen wrote:

   ...

> When I started out in radar signal processing this was routine, and > where I worked it was called a DDA (Digital Differential Analyzer) > instead of a DDS.
A DDA is a more general circuit. A pair can be used to generate oblique lines, circles, and ellipses on an NC milling machine. Ans lots more! Jerry -- Engineering is the art of making what you want from things you can get. &macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;
On Sep 23, 2:52 am, beno...@gmail.com wrote:
> Hi! > > I'm trying to implement a digital chirp signal in hardware. > The formula is: > > c[n] = exp(j*ang[n]), n=0,1,2,...,N-1 > ang[n] = (pi*n^2) / N. > > N is any *prime* number such that N < 1320. > > I have an idea to do that: sampling the exp(j*k) ant insert the > samples into a look up table, > I'll need to calculate an ang[n] every time so c[n] = LUT{ang[n]}. > (LUT = Look UP Table). > > But I don't like this idea because it's very generic and not > efficient. > > I thought about using "chirp z-transform" to generate the digital > signal chirp, but > I don't how... > > Any suggestions? Any other ideas to generate the signal? > > Thanks in advance
Hello, try looking up how the ghost cancellation reference is defined for analog TV. It is a nice bandlimited frequency chirp. It will be easy enough to scale to your application. Clay p.s. google for +tv +gcr you should get a lot of hits.
robert bristow-johnson wrote:
> On Sep 24, 4:46 am, John Monro <johnmo...@optusnet.removethis.com.au> > wrote: >> robert bristow-johnson wrote: >>> On Sep 23, 5:58 pm, John Monro <johnmo...@optusnet.removethis.com.au> >>> wrote: >>>> Hi Benome, >>>> What you are proposing can be quite efficient and is similar to the >>>> standard technique used for Direct Digital Synthesis, except that you >>>> need a progressibely increasing phase-increment instead of a constant >>>> phase-increment. >>> but that's the "LUT" method and he said he didn't like it. (was >>> interpolation the reason the OP didn't like it?) >> Dunno! The Op was looking for suggestions. The fact that I reduced the >> size of the LUT by a factor of eight may make him reconsider using a LUT. >> >> I don't quite understand your references to interpolation, as the method >> I suggested does not use interpolation. >> >> >> >> >> >>>> To generate the phase-increment, you could generate successive n^2 >>>> values and then multiply by each by pi/N as you suggest, but in fact >>>> there is a simpler way which uses second differences and requires no >>>> multiplies. Actually, you don't need to generate n^2 explicitly at >>>> all, as you will see in the following routine:- >>>> CHIRP ROUTINE >>>> Initially: ang =0, D1 = -pi/n, D2 = 2*pi/N >>>> These variables are all fixed-point, with an integer part and a >>>> fractional part. >>>> 1. Look up cos(ang) and sin(ang)to give you the required exp(j*ang). >>>> 2. D1 = D1 + D2 >>>> 3. ang = ang + D1 >>>> 4. If ang < ang(max), go to (1) >>>> 5. End of chirp. >>>> The angle (ang) needs to be held in a long register. The high-order >>>> bits are used to access the LUT and the low-order bits are there to >>>> minimise the accumulated phase error. >>>> By the way, you can reduce the size of the LUT by a factor of eight by >>>> taking advantage of the symmetry of the cos and sin functions. Have a >>>> LUT for only cos(ang), and then only in the range to 0 to pi/2. You >>>> then use appropriate rules to find cos(arg) and sin(arg) for any value >>>> of (arg). >>> one thing to think about is that the chirp signal is not bandlimited >>> (even if the range of the instantaneous frequency of the chirp always >>> had magnitude less than Nyquist). sampling it (which is what you're >>> doing if your interpolation gave you the precise time-domain values of >>> the sin(t^2) function. theoretically just outputting those LUT >>> values, even if the bandlimited interpolation was perfect (using an >>> arbitrarily long sinc(.) function), it's not the same as sampling a >>> "properly" bandlimited chirp. you might have to think about images >>> and aliases if the chirp sweep rate is fast enough. >> I agree that this chirp signal is not bandlimited, but do not understand >> the reference to "...your interpolation.." > > in step #1: > >> 1. Look up cos(ang) and sin(ang)to give you the required exp(j*ang). > > how do you "look up" when ang takes on a value that, when scaled to be > an index into the LUT, is not precisely an integer value (so a single > look-up isn't precise)? what do you do with "the low-order bits > [that] are there to minimise the accumulated phase error"?
The chirp-signal generator algorithm that I described is a minor variation on the usual Direct Digital Synthesis method in that the value of the phase-increment is steadily increasing, rather than staying constant. In DDS, the number of high-order bits is designed to keep the amplitude-noise and phase-noise down to acceptable levels. These bits contribute to the value of the cos or sin functions being generated. The remaining low-order bits in the phase-accumulator register are there to minimise phase and frequency error and don't contribute immediately to the cos or sin function being generated. The low-order bits are not ignored in any way: they influence what gets carried up into the high-order bits when subsequent additions are made to the register. Keeping these low-order bits ensures that the long-term phase and frequency accuracy are a lot better than the short-term phase and frequency accuracy.
> i presume > that step #1 has some interpolation intrinsic to the operation to > handle such values of "ang", since to reference to this is made in the > other steps. now, that's okay John, i'm not criticizing. i'm only > pointing out that whatever you're doing with these "low-order bits", > even simply ignoring them and throwing them away (regarding table > lookup, i know that you keep them in the phase accumulator so that you > have some decent precision on the instantaneous frequency), whatever > reasonable thing you do with those bits in the table lookup (incl. > ignoring them) is a form of interpolation.
OK Robert, I see what you are driving at. I agree that some sort of interpolation could be a good way to reduce the size of the LUT or to increase the accuracy of the cos(ang) and sin(ang) values, but actually I was not assuming that this was taking place. In fact it did not occur to me to mention that possibility. If you did interpolate in the way you describe below I expect that you would still want to maintain a reasonable number of low bits to improve the longer-term phase and frequency accuracy. I realise that in the type of application that the chirp would be used for, the long-term phase and frequency accuracy are less of a concern than they are in the usual DDS application. now let's suppose, for the
> point of illustration is that what you do with those bits is apply > them (and the neighboring values in the table) to an idealized > reconstruction formula like: > > +M/2 > x(m+a) = SUM{ x[(m+n)modulo_N] * sinc(a-n) } > n=-M/2+1 > > where x[n] are the discrete samples of the sine in the table and for > an integer, n, x(n) = x[n] (with no meaningful loss of generality, i'm > normalizing T for simplicity) and M is a very large positive integer > and is the number of terms of the summation in the interpoation. x(t) > is the continuous-time output. m+a is the "ang" parameter where m is > the integer part (m = floor(m+a)) or the "high-order bits" and "a" is > the fractional part (a = (m+a) - ang(m+a) and 0 <= a < 1). also x[n] > are the samples for a single cycle of a sinusoid function: > > x[n] = cos(2*pi*n/N); > > or, conceptually, > > x[n] = e^(j*2*pi*n/N); (j^2 = -1) > > where N is the size of the lookup table and > > (m+n)modulo_N = (m+n) - N*floor((m+n)/N) . > > suppose you do that for interpolation in the LUT. now that original > sinusoid represented in the table is clearly bandlimited to much less > than Nyquist, so the ideal sinc-based reconstruction should give you > pretty much exactly the sinusoid an *any* phase (besides phases that > are 2*pi*n/N) and it shouldn't matter at all what the phase is or how > it is increasing (like a t^2 function). what should come out of > > +M/2 > x(m+a) = SUM{ x[(m+n)modulo_N] * sinc(a-n) } > n=-M/2+1 > > should be awful damn close to > > x(m+a) = e^(j*2*pi*(m+a)/N) > > and, in fact, can come out exact if M=N and we change the sinc(.) to > the Dirichlet thingie: > > +N/2 > x(m+a) = SUM{x[(m+n)mod_N] * sin(pi*(a-n))/(N*sin(pi/N*(a-n)))} > n=-N/2+1 > > (i hope Google Groups didn't wrap that nastily.) > > i'm still debating with myself if, for even N, it should be sin(pi*t)/ > (N*sin(pi/N*t)) or sin(pi*t)/(N*tan(pi/N*t)) for exact > reconconstruction (and i wish someone can settle this and persuade me > of the answer). now, if *that* idealized interpolation is what you > are doing for step #1, then we know we are taking this ideal band- > UNlimited chirp and sampling it at equal spaced times. but since it's > not bandlimited, what are the consequences of sampling it? if your > interpolation is less ideal, it's even a more difficult question to > how nasty the chirp output will be. maybe not very nasty if N is > large and the sweep rate is sufficiently slow. > > r b-j >
As mentioned, there was no interpolation, so the difficulties of interpolating a chirp do not come up. In responding to the OP's query about generating the chirp I did not cover systems issues like band-limiting. Regards, John