Forums

Newbie with questions about decimation for SDR

Started by tomb18 3 years ago27 replieslatest reply 3 years ago1059 views

Hi,

I currently have software that will plot a spectrum from I and Q data obtained from a software defined radio.  The slowest sampling rate is 2msps.  This is fine, I am collecting packets from the device and sending them to an FFT and then eventually plotting them.  I use C# and Mitov's SignalLab.

The thing is that 2msps is too fast and I want to downsample by a factor of 8.  I understand that I need first to use an FIR filter to filter the signal, selecting a cut off frequency and other parameters.  I then just drop 7 samples out of 8.  The question is that everywhere I look it just says to choose the filter parameters but doesn't give any advice on how to do it.

So I think that I need a cutoff frequency of twice the Nyquist which in this case would be 125 kHz.  What about all the other parameters?  

I also heard that I should do this in stages.  For example I could decimate by 2 then 2 the 2 and finally by 2 again.  But, where do I find the algorithms to do so?  Are there example code fragments available for FIR filters?  What about decimation?

Finally, I assume I need to run this code for the I and Q data separately?

Thanks for any help!

#SDR #Resampling

[ - ]
Reply by sgrassiMay 16, 2016

To obtain the filter coefficients without effort, I would use the resample function of Matlab or octave. You specify the resampling factor (p/q) and give the input signal, and it gives you back the resampled signal and the coefficients of the filter used.

www.mathworks.com/help/signal/ref/resample.html

http://octave.sourceforge.net/signal/function/resa...

If you want to know more about resampling I find it very nicely explained in (one of the reference given in octave function)

Ref [1] J. G. Proakis and D. G. Manolakis, Digital Signal Processing: Principles, Algorithms, and Applications, 4th ed., Prentice Hall, 2007. Chap. 6

Hope it helps, 

regards, 

Sara

[ - ]
Reply by tomb18May 16, 2016

Thanks for all the answers. I Guess, there is no absolute rule here and I have to play.  So I started with 32 taps (Mitov has an FIR designer) without downsampling.  The resulting waveform did not have a sharp enough fall off. And the baseline of the spectra was no longer flat over the region of interest. I then went to 64 taps and it looks much better over the frequency range I wish.

So now the next question will really show how much of a newbie I am.

So currently I am taking 16384 samples at a time and feeding them to a FFT with an order of 16384.  This results in a spectrum with a resolution of 2048000/16384 Hz.  Now, when I decimate, I guess I collect 16384 samples from the above (after dropping 7 out of 8 of them) and then feed them to the FFT.  That will give me 2048000 / 8 = 256000 sps and a resolution of 256000/16384 = 15 Hz. So now the spectrum will now show 256000 hz bandwidth with the bins 15 Hz apart. Is this all correct?  

[ - ]
Reply by JOSMay 16, 2016

Also be aware of the sndfile-resample command-line program in the libresample package.  I use that all the time.  It is free open-source so you can study the software implementation. That and other resources are listed at https://ccrma.stanford.edu/~jos/resample/Free_Resa...


[ - ]
Reply by drmikeMay 16, 2016

You have the right idea.  The advantage of doing decimation in stages is that the first few stages can be "sloppy" - you need a few terms with a slow roll off.  Your final stage will need more taps so the cut off is sharper.  But you'll be closer to where you need to be, so it won't take as many taps.  

The other parameters are how good a cut off do you need?  Do you need to be 40 dB down or 120?  That drastically affects the number of taps.  A great place to play around with it is here

You can filter the stream before you break out I and Q, or you can do both separately if that is what you are starting with.  If you do it before it takes half the processing power - and you have decimated the raw data so the rest of your processing will require fewer MIPs.  If you can do it before, it will make your life easier all around.

There are several tricks to decimation - you need to make sure you don't overlap or cutoff information you want.  I think reading a few books will help a lot.  "Understanding Digital Signal Processing" is a good place to start (chapter 7).


[ - ]
Reply by tomb18May 16, 2016

Thanks for all the answers. I Guess, there is no absolute rule here and I have to play.  So I started with 32 taps (Mitov has an FIR designer) without downsampling.  The resulting waveform did not have a sharp enough fall off. And the baseline of the spectra was no longer flat over the region of interest. I then went to 64 taps and it looks much better over the frequency range I wish.

So now the next question will really show how much of a newbie I am.

So currently I am taking 16384 samples at a time and feeding them to a FFT with an order of 16384.  This results in a spectrum with a resolution of 2048000/16384 Hz.  Now, when I decimate, I guess I collect 16384 samples from the above (after dropping 7 out of 8 of them) and then feed them to the FFT.  That will give me 2048000 / 8 = 256000 sps and a resolution of 256000/16384 = 15 Hz. So now the spectrum will now show 256000 hz bandwidth with the bins 15 Hz apart. Is this all correct?  

[ - ]
Reply by Tim WescottMay 16, 2016

All good advice so far, to which I might add a couple of things.

First, it's a good idea to keep in mind why you're filtering (to beat down stuff you don't care about that might alias into your "real" signal) and the anticipated signal content.  This helps in intelligent filter design.  This paper isn't really written with your question in mind (at least your question doesn't suggest a misunderstanding of the Nyquist rate), but you may find it useful even so: http://wescottdesign.com/articles/Sampling/samplin....  There's a section in there specifically on anti-alias filtering which may or may not prove useful.

Second, do a search on polyphase filtering.  The basic idea is that you don't have to run the whole thing through a FIR filter and then discard 7/8 of your work -- instead, you can run the FIR filter for just the output points you're going to keep.  You'll still end up doing a ton of processing, but one ton is better than eight, right?

[ - ]
Reply by tomb18May 16, 2016

Thanks for all the answers. I Guess, there is no absolute rule here and I have to play.  So I started with 32 taps (Mitov has an FIR designer) without downsampling.  The resulting waveform did not have a sharp enough fall off. And the baseline of the spectra was no longer flat over the region of interest. I then went to 64 taps and it looks much better over the frequency range I wish.

So now the next question will really show how much of a newbie I am.

So currently I am taking 16384 samples at a time and feeding them to a FFT with an order of 16384.  This results in a spectrum with a resolution of 2048000/16384 Hz.  Now, when I decimate, I guess I collect 16384 samples from the above (after dropping 7 out of 8 of them) and then feed them to the FFT.  That will give me 2048000 / 8 = 256000 sps and a resolution of 256000/16384 = 15 Hz. So now the spectrum will now show 256000 hz bandwidth with the bins 15 Hz apart. Is this all correct?  

Thanks

[ - ]
Reply by Tim WescottMay 16, 2016

"So now the spectrum will now show 256000 hz bandwidth with the bins 15 Hz apart. Is this all correct?"

Yup.

[ - ]
Reply by tomb18May 16, 2016

Hi,

Where do I reply here? Individually?  Or to my initial post?


[ - ]
Reply by Tim WescottMay 16, 2016

In general, if I'm replying to everyone I reply to my original post (or to the head of a whole sub-thread, if that's appropriate).  I'f I'm going to reply to just post, I put it there.

Use your judgement.  Ask yourself what would make it easiest to read if you were the reader rather than the author, and do that.

[ - ]
Reply by Rick LyonsMay 16, 2016

A simple, and efficient, solution to decimation by 8 is to perform three separate stages of decimation by 2. And in each decimation by 2 stage use a half-band FIR filter. Your best bet is to read the blog at:

https://www.dsprelated.com/showarticle/903.php

You apply the output of the 1st filter to the input of your 2nd filter.  Next you apply the output of the 2nd filter to the input of your 3rd filter. Your final output is the output of your 3rd filter. Remember, each filter operates as follows: Apply two consecutive samples to a filter’s input and compute a single filter output. That way you do not have to perform any sort of decimation operation.  Good Luck.

[ - ]
Reply by napiermMay 16, 2016

For a fixed down-sample by 8 I would also suggest half-band filters.  The nice thing about them is that each one is processing at 1/2 the rate of the sample input.  Also, the 1st one has to run the fastest is also the simplest one.  Note that I'm a hardware guy.  The best solution for SW could be different.

In GNU-Radio there are many down-samplers.  One surprising result there is that the FFT-based FIR filters are faster than the direct implementation for anything over 10-12 taps.  The code is worth looking at, many man-years of work invested.  One of them lets you select which Nyquist zone you want to keep, may not be the one already at base-band.




[ - ]
Reply by tomb18May 16, 2016

Ok, I have it all working....sort of.

For starters I just wanted to get things going simply.  So I am using an FIR filter with 64 taps and a cut off frequency of 125000.

I then get rid of 7 of the 8 samples.

My spectrum is now downsampled and the bandwidth is just as expected and the resolution is about 15 Hz.

However it looks weird.  Here's a picture now:decimated.jpg and here is before (zoomed in)NB-decimation.jpg

In both cases everything is the same.  I use a blackman window.

Any ideas why this seems so funny?

Thanks


[ - ]
Reply by Rick LyonsMay 17, 2016

tomb18, don’t worry. We’ll figure this all out.

First, do not multiply anything by a window sequence, such as a Blackman window. We’ll worry about windowing after we ensure that your filtering/decimation is working.

[1] So that we have an idea of the nature of the signal you’re working with, show us the spectrum of your original signal (Fs = 2.048 Mhz) over the frequency range of 0 Hz –to- 2.048/2 Mhz.

[2] Next, show us the spectrum of your filter output signal (Fs = 2.048 Mhz) over the frequency range of 0 Hz –to- 2.048/2 Mhz.

[3] Finally, show us the spectrum of your decimated filter output signal (Fs = 256 kHz) over the frequency range of 0 Hz –to- 256/2 kHz.

Make sure all your spectral plots extend from zero Hz –to- half the time-domain signal’s Fs sample rate.

[ - ]
Reply by tomb18May 17, 2016

Hi,

Ok, Thanks.

[1] Here is the original showing -2.048MHz to +2.048MHz (it's a quadrature signal) 2msps.jpg

[2] Here is the 2msps with the FIR filter 2msps-with-filter.jpg

[3] Here is the downsampled result without the filter DS-no-filter.jpg

[4] Here is downsampled and filtered. filtered-and-decimated.jpg

I can't provide wider frequency ranges for 3 and 4 since the software will automatically show only the width of the fft.  


[ - ]
Reply by Rick LyonsMay 17, 2016

Ah ha tomb18. You’re making things exciting for us by NOT labeling the horizontal axes of your spectra in Hz. Now we have to guess, “Where is zero Hz on the horizontal axes?” Looking at your first spectrum I decided to assume that zero Hz was at the left end of the horizontal axis. Then, looking at your second spectrum I see that my assumption appears to be wrong! The second spectrum implies that zero Hz is in the center of the horizontal axis. So now I’m confused.

tomb18, delete those four spectral plots and generate new spectral plots with the horizontal axes labeled in Hz. And make sure you know EXACTLY what your FFT and ‘spectral plotting’ software are doing with respect to the horizontal axes’ frequency values in your spectral plots. Make sure you know where zero Hz is located on those axes.

[ - ]
Reply by tomb18May 17, 2016

Hi,

It's up to me to label the axis.,  The software will not do this automatically and I didn't bother at this time since it will depend on the tuning of the radio.  However, 0 Hz is in the middle in all cases.

The output of the FFT is plotted by taking the 0-8191 samples and moving them up to 8192 to 16383 and taking the top samples and moving them down.  The upper two have a span of 2.048MHz and the lower two a span of 256000 Hz.  Just as expected.  All peaks correspond perfectly to the AM broadcast band that I am using.

However, the issue I see is that the downsampled spectra now have huge noise.  And the peaks are no longer look like they should. They look like upside down V's with very slow rise times.  What is wrong there?

[ - ]
Reply by Rick LyonsMay 18, 2016

tomb18, looking at your ‘2msps-with-filter.jpg’ spectrum, if you’re decimating by 8 it appears to me that your lowpass filter’s passband is too wide. For decimation by 8 the center of your filter’s positive-frequency transition region should be located at 128 kHz. Can you prove to yourself that your filter’s positive-frequency transition region is centered at 128 kHz? Your lowpass filter’s positive-frequency passband should extend from zero Hz -to- somewhat less than 128 kHz. Do those words describe your lowpass filter?

[ - ]
Reply by tomb18May 18, 2016

Hi, Thanks.

Here is a picture of the filters response.  filter-response.jpg

Should this be set at 256000 or 128kHz? 

I tried both and it makes no difference to the outcome. The actual zoomed in spectra looks very weird too. zoomed-decimated.jpg

I use the same FFT routine.


[ - ]
Reply by Rick LyonsMay 18, 2016

tomb18, the positive-frequency cutoff frequency of your lowpass filter should be just less than +128 kHz.

If you pass your original input signal through your filter, the output freq range of –128 kHz to +128 kHz spectrum of your filter should be almost identical to the freq range of –128 kHz to +128 kHz spectrum of the signal applied to your filter’s input.Are you seeing that?

Your spectrum in ‘zoomed-decimated.jpg’ has no meaning for me because you have not labeled the horizontal axis in Hz. tomb18, I’ll be gentle. Please do not expect a plot of a mysterious undefined red squiggly curve to have meaning for us.

[ - ]
Reply by tomb18May 18, 2016

Hi,

Yes I am seeing that the -128 KHz to +128 kHz is pretty well identical to the filter output.

Where I see the issue is that the 256kHz downsampled spectrum looks very different than the 256kHz range of the original spectrum.

Here is the 2msps spectrum zoomed in to a 200 kHz area: zoomed-in-2msps.jpg

Here is the same region after dopwnsampling: downsampled.jpg

I guess the question is why does it look like this? It shouldn't. Why is the noise between signals so large?

As to the above "zoomed-decimated.jpg" lack of meaning I guess I should have explained this mode.  Each point is approximately 15Hz apart from the others. I did this to show that there appears to be a weird pattern.  Is this why the decimated spectrum looks so funny?

As far as I can tell, everything I am doing is OK.  But the final spectrum looks wrong.

[ - ]
Reply by Rick LyonsMay 19, 2016

Hi tomb18.  I said "we'll figure this out", and I meant it.  Is there a way I can download your original 2.048-MHz sample rate file in the 'xxx.wav' format? (That way I can use MATLAB to experiment with your original signal.)   Can you e-mail such a file to me, or perhaps place it in a "Dropbox" so I can download it?

My e-mail address is: R.Lyons@ieee.org

[After we solve this problem I'd sure like to know what hardware you used to capture your AM radio signal.  I bought one of those el-cheapo RTL-SDR receivers but it won't tune down to the AM broadcast band as I would like.]

[ - ]
Reply by tomb18May 19, 2016

Hi,

I'm using the SDRPlay device. It has a range of 100 kHz to 2 gHz with a bandwidth of 8 MHz.  One of the least expensive "high end" SDRs.

As to a WAV file, I can't.  I'm starting from actually reading the packets from the device.  I suppose I could write some code to do so but that will take a while.


[ - ]
Reply by tomb18May 19, 2016

Problem is solved...

The SDRPlay device provides packets with different numbers of samples depending on the sampling rate.  Unfortunately, these packets are not based on a power of 2.  So at 2msps, the packets have 332 samples in them. For my tools, I need a packet size of 16384 for the FFT.  So I need 48 packets with a remainder of 256 samples to make 16384 samples. So this is what I was doing.

However, there were 80 samples left over and therefore missing from the next 16384 sample block, (The "leftovers") that I was not using in error.  When I sampled at 2msps, this had no effect.  However when downsampling the missing 80 samples out of 256000 made such a difference that is "smeared" the results!

None the less, this has been an excellent learning experience for me.  I have never downsampled before and running through the steps taught me a lot.  Now to concentrate on the half band FIR filter.

Many thanks for your help.  If you want more details on the SDRPlay let me know with email at tomb18 at videotron dot com.


[ - ]
Reply by Rick LyonsMay 19, 2016

tomb18, thanks for the pointer to the SDRPlay hardware!!

[ - ]
Reply by jyaMay 16, 2016

In the spirit of "First make it work, then make it pretty" I'll point out that you needn't actually calculate the samples that you intend to discard.

Jerry

[ - ]
Reply by tomb18May 16, 2016

Yes, I got that.  But why do the spectra look so different?