all4dsp <all4dsp@comcast.net> wrote: (snip)> I'm not sure if FFT is best either. I also thought about using simple > polynomial interpolation or other type of interpolation instead of > sinc-based (more accurate than linear interpolation) to find zero crossings > directly from the ADC array, but there are so many ways to interpolate a > curve through a few points near the zero-crossing, I don't have a good > enough grasp of the theory to guide me as to which method would work in > this situation. At least the sinc interpolation relies on the frequency > content of the entire waveform to decide where the new points should be > located, so I see it as the gold standard for this sort of thing (although, > unfortunately I don't need 99% of the interpolated points this method > returns).The sinc interpolation does use the frequency content, but then you do linear interpolation, which probably loses some of that. Many people would do cubic spline interpolation. If you need to do better, then a higher odd order polynomial could be used. The wikipedia entry for Cubic spline has a good explanation. There is also Numerical Recipes for explanation and source code. -- glen
conventional wisdom how to upsample very large arrays, accurately?
Started by ●January 6, 2010
Reply by ●January 6, 20102010-01-06
Reply by ●January 6, 20102010-01-06
On Wed, 06 Jan 2010 06:36:03 -0600, "all4dsp" <all4dsp@comcast.net> wrote:>Hello all, > >I'd like to know what methods are available, if any, for upsampling a >large array. I have an application where I need to sinc interpolate an >array of 1E9 (billion) double precision data points output by an ADC by 4x. >The array captures a clock waveform centered at zero volts and is >oversampled (meets Nyquist). Don't know if it matters, but my goal is >simply to very accurately find the zero-point crossings. By upsampling a >factor of 4 times, I'll then use linear interpolation to get the final >number. The point is, I really only need very accurate data during the >rise/fall times (not sure if this opens a new degree of freedom). Sample >spacing is periodic. I can throw out the initial and final 5% or so of the >array if needed to eliminate artifacts, but I need the majority of data in >the middle. > >My current system takes 10 hours to computer a 1E9 FFT (actually, size is >2^30, or 1.074G points). The CPU only works ~1%, while virtual memory is >57GB. Not sure I could stuff 3G pts into the FFT then IFFT it. Would like >to get it down to as short as possible (10 minutes?). > >I see the Overlap-Add method in Chapter 18 of dspguide.com handles large >arrays, but I'm not sure how to use it for upsampling. > >Rick Lyons excellent blog on Computing Large DFTs Using Small FFTs >(http://www.dsprelated.com/showarticle/63.php) shows how to get an DFT >using smaller FFTs, but I'm not sure how to do the upsampling. > >Could someone help me understand what the conventional wisdom is to do >this? I'm not a DSP expert, although I'm a fast study, if you could >reference literature I could access (no research papers where someone >created an algorithm no one uses though) I'd be very much appreciative. >Thanks in advance.Hello all4dsp, If your input signal (the signal that's to be interpolated) is a "squarewave-like" signal, then maybe the linear interpolation method in: [1] R. Losada and R. lyons, "Reducing CIC Filter Complexity", IEEE Signal processing magazine, July 2006, pp.124-126 would be of some help to you. That article describes a linear time-domain interpolation scheme that requires no multiplications regardless of the interpolation factor. A slightly expanded version of the material in the above article can be found in Chapter 6 in: [2] R. Lyons, Editor, "Streamlining Digital Signal Processing, A Tricks of the Trade Guidebook", IEEE-Wiley Interscience Publishing, 2007. all4dsp, if you're unable to gain access to the above article or the above book chapter, then let me know. Good Luck, [-Rick-]
Reply by ●January 6, 20102010-01-06
On Wed, 06 Jan 2010 06:36:03 -0600, "all4dsp" <all4dsp@comcast.net> wrote: [Snipped by Lyons]> >I see the Overlap-Add method in Chapter 18 of dspguide.com handles large >arrays, but I'm not sure how to use it for upsampling. > >Rick Lyons excellent blog on Computing Large DFTs Using Small FFTs >(http://www.dsprelated.com/showarticle/63.php) shows how to get an DFT >using smaller FFTs, but I'm not sure how to do the upsampling.Hi all4dsp, If your input data block length (and FFT size) is N, then you can interpolate by four by: (1) Perform your initial N-point FFT to compute the X(m) freq samples. (Freq index m goes from 0 -to- N-1.) (2) Create a new freq-domain sequence comprising: the first N/2 X(m) samples, followed by 3*N zero-valued samples, followed by the last N/2 X(m) samples. This step merely inserts 3*N zero-valued samples right in the middle of your original FFT's X(m) samples. (3) Perform a 4N-point inverse FFT on this "expanded" freq-domain sequence to obtain your interpolated-by-four time-domain sequence. I'm assuming that the initial FFT samples, X(m), are *VERY* low in magnitude in the vicinity of Fs/2 Hz (Fs is the original sample rate in Hz). Of course the above scheme's gonna take a super long time to execute because your working with such unbelievably humungous (is that a word?) large data block lengths. Your problem give a whole new meaning to the phrase "number crunching". all4dsp, thinking about this, I'm not at all sure that the above freq-domain interpolation is the most computationally-efficient thing to do. Time-domain interpolation requires digital filtering, and if the number of filter taps is not too large maybe time-domain interpolation is the "way to go". If so, then polyphase filtering methods should be used. all4dsp, ya' know what I'd do if I were you? I'd carefully double-check the need (the necessity) for working with such large-sized blocks of data. And another thing, as Veronica said in the movie "The Fly", I'd "Be afraid. Be very afraid" of numerical errors in performing such large-sized FFTs. In those FFTs you're expecting the software to accuarately estimate the sine and cosine of angles in the range of 360/2^30 degrees!!! Are sine and cosine approximation algorithms (in the FFT software) accurate enough for such super-small angles??? Good Luck, [-Rick-]
Reply by ●January 7, 20102010-01-07
On Jan 6, 4:36�am, "all4dsp" <all4...@comcast.net> wrote:> I'd like to know what methods are available, if any, for upsampling a > large array. I have an application where I need to sinc interpolate an > array of 1E9 (billion) double precision data points output by an ADC by 4x. > The array captures a clock waveform centered at zero volts and is > oversampled (meets Nyquist). Don't know if it matters, but my goal is > simply to very accurately find the zero-point crossings. By upsampling a > factor of 4 times, I'll then use linear interpolation to get the final > number. The point is, I really only need very accurate data during the > rise/fall times (not sure if this opens a new degree of freedom). Sample > spacing is periodic. I can throw out the initial and final 5% or so of the > array if needed to eliminate artifacts, but I need the majority of data in > the middle. > > My current system takes 10 hours to computer a 1E9 FFT ...What is the density of zero crossings to samples? If it's low enough, you can do an initial linear or quadratic interpolation, then a sinc interpolation (polyphase, etc.) only surrounding these first estimate points. Interpolate again to get a second estimate (and possibly rinse-and-repeat until you converge sufficiently, if needed). YMMV. -- rhn A.T nicholson d.0.t C-o-M http://www.nicholson.com/dsp.html
Reply by ●January 7, 20102010-01-07
On 6 Jan., 18:40, Clay <c...@claysturner.com> wrote:> On Jan 6, 7:36�am, "all4dsp" <all4...@comcast.net> wrote: > > > > > > > Hello all, > > > I'd like to know what methods are available, if any, for upsampling a > > large array. I have an application where I need to sinc interpolate an > > array of 1E9 (billion) double precision data points output by an ADC by 4x. > > The array captures a clock waveform centered at zero volts and is > > oversampled (meets Nyquist). Don't know if it matters, but my goal is > > simply to very accurately find the zero-point crossings. By upsampling a > > factor of 4 times, I'll then use linear interpolation to get the final > > number. The point is, I really only need very accurate data during the > > rise/fall times (not sure if this opens a new degree of freedom). Sample > > spacing is periodic. I can throw out the initial and final 5% or so of the > > array if needed to eliminate artifacts, but I need the majority of data in > > the middle. > > > My current system takes 10 hours to computer a 1E9 FFT (actually, size is > > 2^30, or 1.074G points). The CPU only works ~1%, while virtual memory is > > 57GB. Not sure I could stuff 3G pts into the FFT then IFFT it. Would like > > to get it down to as short as possible (10 minutes?). > > > I see the Overlap-Add method in Chapter 18 of dspguide.com handles large > > arrays, but I'm not sure how to use it for upsampling. > > > Rick Lyons excellent blog on Computing Large DFTs Using Small FFTs > > (http://www.dsprelated.com/showarticle/63.php) shows how to get an DFT > > using smaller FFTs, but I'm not sure how to do the upsampling. > > > Could someone help me understand what the conventional wisdom is to do > > this? I'm not a DSP expert, although I'm a fast study, if you could > > reference literature I could access (no research papers where someone > > created an algorithm no one uses though) I'd be very much appreciative. > > Thanks in advance. > > Here is a way. > > http://www.claysturner.com/dsp/DSPinterpolation.pdf > > This can be coded up and runs reasonably well. If the "spectral > occupancy" is a fair bit less than Nyquist then the sync kernal can be > replaced by the quicker decaying interpolation function.Hi Clay The problem with your method is to find the zero-crossings in the sampled data. It is well possible that zero-crossings can't be found by multiplying two consecutive numbers and checking the sign of the result. This effect does not depend on the "spectral occupancy". To see this, consider sampling x(t) = sin(pi f t) + sin(pi (1-f) t) for arbitrary f (however, the effect is very nice for f close to 0 or 1) with sampling period T = 1. Your algorithm will not find any zero- crossing. An easy modification to your searching criterion can by circumvented by adding a small offset to x (in this case, only about half of the zero crossings will be found). Regards, Andor
Reply by ●January 7, 20102010-01-07
On Jan 7, 3:31�am, Andor <andor.bari...@gmail.com> wrote:> On 6 Jan., 18:40, Clay <c...@claysturner.com> wrote: > > > > > > > On Jan 6, 7:36�am, "all4dsp" <all4...@comcast.net> wrote: > > > > Hello all, > > > > I'd like to know what methods are available, if any, for upsampling a > > > large array. I have an application where I need to sinc interpolate an > > > array of 1E9 (billion) double precision data points output by an ADC by 4x. > > > The array captures a clock waveform centered at zero volts and is > > > oversampled (meets Nyquist). Don't know if it matters, but my goal is > > > simply to very accurately find the zero-point crossings. By upsampling a > > > factor of 4 times, I'll then use linear interpolation to get the final > > > number. The point is, I really only need very accurate data during the > > > rise/fall times (not sure if this opens a new degree of freedom). Sample > > > spacing is periodic. I can throw out the initial and final 5% or so of the > > > array if needed to eliminate artifacts, but I need the majority of data in > > > the middle. > > > > My current system takes 10 hours to computer a 1E9 FFT (actually, size is > > > 2^30, or 1.074G points). The CPU only works ~1%, while virtual memory is > > > 57GB. Not sure I could stuff 3G pts into the FFT then IFFT it. Would like > > > to get it down to as short as possible (10 minutes?). > > > > I see the Overlap-Add method in Chapter 18 of dspguide.com handles large > > > arrays, but I'm not sure how to use it for upsampling. > > > > Rick Lyons excellent blog on Computing Large DFTs Using Small FFTs > > > (http://www.dsprelated.com/showarticle/63.php) shows how to get an DFT > > > using smaller FFTs, but I'm not sure how to do the upsampling. > > > > Could someone help me understand what the conventional wisdom is to do > > > this? I'm not a DSP expert, although I'm a fast study, if you could > > > reference literature I could access (no research papers where someone > > > created an algorithm no one uses though) I'd be very much appreciative. > > > Thanks in advance. > > > Here is a way. > > >http://www.claysturner.com/dsp/DSPinterpolation.pdf > > > This can be coded up and runs reasonably well. If the "spectral > > occupancy" is a fair bit less than Nyquist then the sync kernal can be > > replaced by the quicker decaying interpolation function. > > Hi Clay > > The problem with your method is to find the zero-crossings in the > sampled data. It is well possible that zero-crossings can't be found > by multiplying two consecutive numbers and checking the sign of the > result. This effect does not depend on the "spectral occupancy". To > see this, consider sampling > > x(t) = sin(pi f t) + sin(pi (1-f) t) > > for arbitrary f (however, the effect is very nice for f close to 0 or > 1) with sampling period T = 1. Your algorithm will not find any zero- > crossing. An easy modification to your searching criterion can by > circumvented by adding a small offset to x (in this case, only about > half of the zero crossings will be found). > > Regards, > Andor- Hide quoted text - > > - Show quoted text -Certainly true, I did this for an application years ago where the received data had been substantially band passed. Zero crossing location is used for ranging lightning. If I had to do this today, I would use better methods involving optimized interpolation kernals. Thanks for the input, Clay
Reply by ●January 7, 20102010-01-07
Thanks everyone for your replies, you've given me a lot to think about. There are obviously a lot of tradeoffs in selecting a method.>I'm assuming that the initial FFT samples, X(m), are *VERY* >low in magnitude in the vicinity of Fs/2 Hz (Fs is the >original sample rate in Hz).This is true.>In those FFTs you're expecting the software to accuarately >estimate the sine and cosine of angles in the range >of 360/2^30 degreesWow, good point! It sounds like polyphase filtering might be the best choice.
Reply by ●January 7, 20102010-01-07
On 7 Jan., 16:42, Clay <c...@claysturner.com> wrote:> On Jan 7, 3:31�am, Andor <andor.bari...@gmail.com> wrote: > > > > > > > On 6 Jan., 18:40, Clay <c...@claysturner.com> wrote: > > > > On Jan 6, 7:36�am, "all4dsp" <all4...@comcast.net> wrote: > > > > > Hello all, > > > > > I'd like to know what methods are available, if any, for upsampling a > > > > large array. I have an application where I need to sinc interpolate an > > > > array of 1E9 (billion) double precision data points output by an ADC by 4x. > > > > The array captures a clock waveform centered at zero volts and is > > > > oversampled (meets Nyquist). Don't know if it matters, but my goal is > > > > simply to very accurately find the zero-point crossings. By upsampling a > > > > factor of 4 times, I'll then use linear interpolation to get the final > > > > number. The point is, I really only need very accurate data during the > > > > rise/fall times (not sure if this opens a new degree of freedom). Sample > > > > spacing is periodic. I can throw out the initial and final 5% or so of the > > > > array if needed to eliminate artifacts, but I need the majority of data in > > > > the middle. > > > > > My current system takes 10 hours to computer a 1E9 FFT (actually, size is > > > > 2^30, or 1.074G points). The CPU only works ~1%, while virtual memory is > > > > 57GB. Not sure I could stuff 3G pts into the FFT then IFFT it. Would like > > > > to get it down to as short as possible (10 minutes?). > > > > > I see the Overlap-Add method in Chapter 18 of dspguide.com handles large > > > > arrays, but I'm not sure how to use it for upsampling. > > > > > Rick Lyons excellent blog on Computing Large DFTs Using Small FFTs > > > > (http://www.dsprelated.com/showarticle/63.php) shows how to get an DFT > > > > using smaller FFTs, but I'm not sure how to do the upsampling. > > > > > Could someone help me understand what the conventional wisdom is to do > > > > this? I'm not a DSP expert, although I'm a fast study, if you could > > > > reference literature I could access (no research papers where someone > > > > created an algorithm no one uses though) I'd be very much appreciative. > > > > Thanks in advance. > > > > Here is a way. > > > >http://www.claysturner.com/dsp/DSPinterpolation.pdf > > > > This can be coded up and runs reasonably well. If the "spectral > > > occupancy" is a fair bit less than Nyquist then the sync kernal can be > > > replaced by the quicker decaying interpolation function. > > > Hi Clay > > > The problem with your method is to find the zero-crossings in the > > sampled data. It is well possible that zero-crossings can't be found > > by multiplying two consecutive numbers and checking the sign of the > > result. This effect does not depend on the "spectral occupancy". To > > see this, consider sampling > > > x(t) = sin(pi f t) + sin(pi (1-f) t) > > > for arbitrary f (however, the effect is very nice for f close to 0 or > > 1) with sampling period T = 1. Your algorithm will not find any zero- > > crossing. An easy modification to your searching criterion can by > > circumvented by adding a small offset to x (in this case, only about > > half of the zero crossings will be found). > > > Regards, > > Andor- Hide quoted text - > > > - Show quoted text - > > Certainly true, I did this for an application years ago where the > received data had been substantially band passed. Zero crossing > location is used for ranging lightning. If I had to do this today, I > would use better methods involving optimized interpolation kernals.I don't think the interpolation kernels are the problem (although they certainly can be). I think the main problem is to find the rough location of the zero- crossing in the sampled data. As my example shows, a zero-crossing of the underlying sampled function will not necessarily show up as a "zero-crossing" in the samples. I would opt for Ron's idea to use quick polynomial interpolation with order >= 3 (can be implemented with very small kernel FIR filters) to oversample the data by 4 or so to find "zero-crossings" in the samples and then use the method of sinc-interpolation outlined in your paper to find the exact location. This reminds me of the old comp.dsp debate on the maximum number of zero-crossings that a bandlimited function may have on any given interval. As it turned out (I think Matt Timmermans came up with the example of the prolate spheroidal wave functions), the number of zero- crossings is not limited ...> > Thanks for the input,Thanks for the paper :-). Regards, Andor
Reply by ●January 7, 20102010-01-07
On 7 Jan., 16:48, "all4dsp" <all4...@comcast.net> wrote:> Thanks everyone for your replies, you've given me a lot to think about. > There are obviously a lot of tradeoffs in selecting a method. > > >I'm assuming that the initial FFT samples, X(m), are *VERY* > >low in magnitude in the vicinity of Fs/2 Hz (Fs is the > >original sample rate in Hz). > > This is true. > > >In those FFTs you're expecting the software to accuarately > >estimate the sine and cosine of angles in the range > >of 360/2^30 degrees > > Wow, good point! It sounds like polyphase filtering might be the best > choice.Since you are upsampling by 4 (an integer factor) you don't even need to implement a polyphase resampler. All you need is a filterbank with three filters. As I wrote in another post, I would suggest a two step procedure: 1. Interpolate by 4 with polynomial FIR filters (very small kernels, e.g. 4 taps for third-order polynomial interpolation): x(n) ---+-----------------> x(n) | | +---------+ |-->| H1(z) |---> x1(n) | +---------+ | | +---------+ |-->| H2(z) |---> x2(n) | +---------+ | | +---------+ +-->| H3(z) |---> x3(n) +---------+ H1, H2 and H3 are the interpolation filters (polynomial filters are simple to design from scratch), your upsampled sequence will be: (... x(n),x1(n),x2(n),x3(n),x(n+1),x1(n+1),x2(n+1),x3(n+1), ... ) Using a third-order polynomial filter requires 3*4*n MACS to resample by 4. For n = 1e9 you have 12e9 MACS, which should take in a couple of seconds on a normal PC. 2. Use a simple method to find zero-crossings in the upsampled data and then "polish" the zero-crossing location with more accurate interpolation - this is outlined here: http://www.claysturner.com/dsp/DSPinterpolation.pdf Regards, Andor
Reply by ●January 8, 20102010-01-08
On 7 Jan., 17:18, Andor <andor.bari...@gmail.com> wrote:> On 7 Jan., 16:48, "all4dsp" <all4...@comcast.net> wrote: > > > Thanks everyone for your replies, you've given me a lot to think about. > > There are obviously a lot of tradeoffs in selecting a method. > > > >I'm assuming that the initial FFT samples, X(m), are *VERY* > > >low in magnitude in the vicinity of Fs/2 Hz (Fs is the > > >original sample rate in Hz). > > > This is true. > > > >In those FFTs you're expecting the software to accuarately > > >estimate the sine and cosine of angles in the range > > >of 360/2^30 degrees > > > Wow, good point! It sounds like polyphase filtering might be the best > > choice. > > Since you are upsampling by 4 (an integer factor) you don't even need > to implement a polyphase resampler. All you need is a filterbank with > three filters. As I wrote in another post, I would suggest a two step > procedure: > > 1. Interpolate by 4 with polynomial FIR filters (very small kernels, > e.g. 4 taps for third-order polynomial interpolation): > > x(n) ---+-----------------> x(n) > � � � � | > � � � � | � +---------+ > � � � � |-->| �H1(z) �|---> x1(n) > � � � � | � +---------+ > � � � � | > � � � � | � +---------+ > � � � � |-->| �H2(z) �|---> x2(n) > � � � � | � +---------+ > � � � � | > � � � � | � +---------+ > � � � � +-->| �H3(z) �|---> x3(n) > � � � � � � +---------+ > > H1, H2 and H3 are the interpolation filters (polynomial filters are > simple to design from scratch), your upsampled sequence will be: > > (... x(n),x1(n),x2(n),x3(n),x(n+1),x1(n+1),x2(n+1),x3(n+1), ... ) > > Using a third-order polynomial filter requires 3*4*n MACS to resample > by 4. For n = 1e9 you have 12e9 MACS, which should take in a couple of > seconds on a normal PC. > > 2. Use a simple method to find zero-crossings in the upsampled data > and then "polish" the zero-crossing location with more accurate > interpolation - this is outlined here:http://www.claysturner.com/dsp/DSPinterpolation.pdfhmm. For fun I implemented this scheme and applied it to the samples of the signal x(t) = sin(pi f t) + sin(pi (1-f) t) + eps s((1-f)/2 t) for 0.5 < f < 1.0, eps > 0 and s(t) the square wave function with period 1, sampled with sampling interval T = 1. The sampling sequence x [n] = x(n) has the property that all samples from the time interval 0 <= t < 1/(1-f) are positive, and all samples from the interval 1/(1-f) <= t < 2/(1-f) are negative. For small eps, the number of zero-crossings of x(t) in both intervals is about 1/(1-f)/f. Unfortunately, upsampling by a factor of 4 with third-order polynomial FIR did not change the fact that the first half of the samples are all positive and the second half are all negative (meaning we can't find the zero-crossings from the sampled signal using the "multiply-two- consecutive-samples-and-check-the-sign" trick). So I think the OP might possibly not get around to using high-order sinc-interpolation filters for the factor 4 upsampling on the whole data set to find the zero-crossings. Regards, Andor
