Re: [Discuss-gnuradio] DDC and Polyphase Channelizer

[Tom Rondeau]

On Wed, Dec 29, 2010 at 4:32 AM, Jimmy Richardson <address@hidden> wrote:
> Hi, Tom:
>
> Please see my comments below.
>
> Thanks
>
> On 12/29/2010 12:58 AM, Tom Rondeau wrote:
>>
>> On Mon, Dec 27, 2010 at 3:48 AM, Jimmy Richardson<address@hidden>
>>  wrote:
>>>
>>> Hi,:
>>>
>>> I'm learning about polyphase channelizer, and I have been reading fred
>>> harris' paper "Digital Receivers and Transmitters Using Polyphase Filter
>>> Banks for Wireless Communications". I have also played with the Matlab
>>> programs at the end of the paper. What I'm trying to do next is to use
>>> GNURadio to do the channelization. I made two attempts, one uses
>>> freq_xlating_fir_filter_ccf to implement a DDC channelizer, another uses
>>> pfb_channelizer_ccf, running against the data in appendix III of harris'
>>> paper. I'm hoping the channelizer result from GNURadio code would match
>>> the
>>> result from Matlab, but I failed on both, I wonder if someone could help
>>> me
>>> spot the problem.
>>>
>>> 1. DDC Channelizer:
>>
>> Using the frequency xlating filters here is not recommended. You'll
>> need to use much longer filters than when using the channelizer. That
>
> Could you clarify? Does this mean I need to change the parameter to
> gr.firdes.low_pass, maybe reduce _cutoff_freq?


I was just saying that, in general, you'll be doing much more work
this way. But since you say later that you're only using this for
comparisons, then you already know that.


>> said, I'm not sure what would cause the problem you identified below.
>> In my opinion, it's not worth doing it this way.
>
> Ok, the reason I did this is because I want a way to verify the result from
> pfb channelizer is correct, since I'm not sure if I'm using it correctly.
>
>>
>>> 1.1 Code: The following is the core part of the code (input parsing and
>>> various other stuff are omitted), the value for the variables can be
>>> found
>>> in 1.2
>>>        self._filter_taps = gr.firdes.low_pass(1, self._input_rate,
>>> self._cutoff_freq, self._trans_width,
>>>
>>> window=gr.firdes.WIN_BLACKMAN_hARRIS)
>>>        print "Number of taps:   _filter_taps = ", len(self._filter_taps)
>>>
>>>        # Blocks
>>>        self.ddc_filter = list()
>>>        # i = 0 to M-1
>>>        for i in xrange(self._M):
>>>            offset = self._channel_bandwidth * (float(self._M)/2 - i)
>>>            print "Channel ", i, "'s offset (Hz): ", offset
>>>
>>>
>>>  
>>> self.ddc_filter.append(gr.freq_xlating_fir_filter_ccf(self._resampling_ratio,
>>>
>>>  self._filter_taps, offset, self._input_rate))
>>>            self.connect(self.head, self.ddc_filter[i],
>>> self.output_files[i])
>>>
>>> 1.2 Variable values:
>>> Input sampling freq (Hz):  _input_rate =  1120000000.0
>>> Resampling ratio:          _resampling_ratio =  28
>>> Output sampling freq (Hz): _output_rate =  40000000.0
>>> Transition width (Hz):     _trans_width =  2000000.0
>>> Number of channels:        _M =  40
>>> Number of samples:         _N =  11200
>>> Oversample rate:           _oversample_rate =  1.42857142857
>>> Channel bandwidth (HZ):    _channel_bandwidth =  28000000.0
>>> Cutoff freq (HZ):          _cutoff_freq =  13000000.0
>>> Number of taps:   _filter_taps =  561
>>> Channel  0 's offset (Hz):  560000000.0
>>> Channel  1 's offset (Hz):  532000000.0
>>> ...
>>> 1.3 Result: The result from GNURadio code is saved to output files, one
>>> for
>>> each channel. Then the results time-domain are plotted in Matlab (-10 to
>>> +10, time 1 to 200), and compared to the result plot from harris' code,
>>> they
>>> are similar, but not the same. Channel 0 to 5, 17, 18, 34 to 39 do not
>>> have
>>> any strong signal, this is reflected in both plots, channel 20's signal
>>> is
>>> exactly the same, but all other channels' plot are different.
>>>
>>> 2. Polyphase Channelizer:
>>> 2.1 Code:
>>>        self._filter_taps = gr.firdes.low_pass(1, self._input_rate,
>>> self._cutoff_freq, self._trans_width,
>>>
>>>  window=gr.firdes.WIN_BLACKMAN_hARRIS)
>>>        # Calculate the number of taps per channel for our own information
>>>        self._taps_per_channel = scipy.ceil(float(len(self._filter_taps))
>>> /
>>>  float(self._M))
>>>        print "Number of taps:   _filter_taps = ", len(self._filter_taps)
>>>        print "Taps per channel: _taps_per_channel = ",
>>> self._taps_per_channel
>>>
>>>        # Construct the channelizer filter
>>>        self.pfb = blks2.pfb_channelizer_ccf(self._M, self._filter_taps,
>>> self._oversample_rate)
>>>
>>>        # Connect the blocks
>>>        self.connect(self.head, self.pfb)
>>>        # i = 0 to M-1
>>>        for i in xrange(self._M):
>>>            self.connect((self.pfb, i), self.output_files[i])
>>>
>>> 2.2 Variable values:
>>> Input sampling freq (Hz):  _input_rate =  1120000000.0
>>> Resampling ratio:          _resampling_ratio =  28
>>> Output sampling freq (Hz): _output_rate =  40000000.0
>>> Transition width (Hz):     _trans_width =  2000000.0
>>> Number of channels:        _M =  40
>>> Number of samples:         _N =  11200
>>> Oversample rate:           _oversample_rate =  1.42857142857
>>> Channel bandwidth (HZ):    _channel_bandwidth =  28000000.0
>>> Cutoff freq (HZ):          _cutoff_freq =  13000000.0
>>> Number of taps:   _filter_taps =  561
>>> Taps per channel: _taps_per_channel =  15.0
>>>
>>> 2.3 Result: The GNURadio program failed with:
>>> "terminate called after throwing an instance of 'std::invalid_argument'
>>>  what():  gr_pfb_channelizer: oversample rate must be N/i for i in [1, N]
>>>
>>> This application has requested the Runtime to terminate it in an unusual
>>> way.
>>> Please contact the application's support team for more information."
>>>
>>> This doesn't seem to make sense, since in this case N = 40, oversample
>>> rate
>>> =1.42857142857, i = 28
>>>
>>> Thanks
>>
>> The "oversample_rate" is not exactly what you are specifying here. The
>> default behavior of this PFB channelizer implementation is to be
>> critically sampled; that is, the sample rate equals the bandwidth. For
>> many applications, you want to have the output sample rate be some
>> other rate, most notably, 2x the bandwidth. The oversample_rate allows
>> you to specify a number real number, but this cannot be an arbitrary
>> number in this case. Because of the way the samples are fed to the
>> filters, the oversample rate, O, has to be set such that the number of
>> filters, N, can be evenly divided by an integer to to get O; that is,
>> N/i = O for i in [1,N]. Take the example in in fred's paper where he
>> has 64 channels and a 1/48 resampling rate. So he's using a 48-to-1
>> downsampling in 64 channels, so his rate is 4/3, which means that the
>> output sampling rate is (4/3)*bw, where bw is the channel bandwidth.
>
> Sorry, false alarm here, it looks like a problem with MinGW's compiler, it
> didn't compile double comparison with 0.0 correctly or something. I tried
> the program under Ubuntu and the error disappeared. I think the number I
> used is the correct number since in my case N=40, i = 28, and O = N/i =
> 1.42857142857

Ok. I didn't calculate the number myself, just assumed it was wrong if
you were being given that error message :)

On a side note, I don't like that problem happening on MinGW.
Honestly, doing an == on floating point numbers isn't a good idea, so
I may need to rethink that check.

>> If you specified that the rate was 2.0, you would feed 32 of the 64
>> channels and have a sample rate that is twice the bandwidth. I notice
>> in the remarks in fred's paper that he's outputting the channels at 2
>> samples/symbol, so the oversample_rate here should be 2.
>
> I'm using fred's example in appendix III, not the 64 channels example. The
> comment summaries the example as below:

Yeah, I know. I was using fred's 64-channel example from the paper
because the numbers are easier to handle/explain.

> % receiver_40z is a demo of a 40 channel receiver, demodulating 30 channels,
> % of nominal symbol rate 20 MHz, separated by 28 MHz centers (1.4 times
> symbol rate)
> % input sample rate is 40*28 = 1120 MHz.
> % receiver performs a 40 point transform on the output of a 40-stage
> polyphase filter
> % the polyphase filter operates at input rate but outputs at 2
> samples/symbol or
> % 40 MHz. The resampling rate is 1120/40 = 28-to-1, thus output from
> 40-channels are
> % computed once for every 28 input samples. channelizer is not matched
> filter,
> % prototype filter is 10% wider than two sided bandwidth of input signal to
> accommodate
> % frequency uncertainty of separate channel centers.
>
>> I hope that makes sense and you can figure out from there what you
>> need to do to get it to work. I've worked with fred for a while now,
>> but I find his Matlab to be almost unreadable (I've told him this
>> before), so I can't quite tell you from a quick look over the code he
>> has in this paper what's really going on.
>
> I got passed the error, thanks, but still have problems with channelizer,
> see below...
>
>> Also, keep in mind that fred always develops he own filters straight
>> from the remez algorithm. If you are using gr.firdes.low_pass, you're
>> not going to get the same output filters. Closer would be to use the
>> blcks2.optfir.low_pass, which uses the PM algorithm (which in turn
>> uses Remez), but I couldn't give you the specs you need to build the
>> same _exact_ filter as in fred's paper. But you'll get close enough
>> that the basic shapes of the signals will be the same, if not every
>> point.
>
> Thanks for the suggestion, I modified filter statement to:
>
>        self._filter_taps = blks2.optfir.low_pass(1, self._input_rate,
>                                                  self._cutoff_freq -
> self._trans_width/2, self._cutoff_freq + self._trans_width/2,
>                                                  0.9, 60, 400)
>
> Note I have to use 400 as extra taps otherwise remez won't work (too many
> extremals -- cannot continue)

Hm, that looks wrong. You shouldn't need to use the extra taps input
just to get the algorithm to converge. It's really meant to be able to
give you minor control over the number of taps you use, specifically
if you need a filter with an even or odd number of taps. Using 400 is
way too many.

The real problem I think comes from your in-band ripple or 0.9. That's
a value in dB and specifies the magnitude of the peak-to-peak
fluctuations in the passband. I was able to use your settings above
but specified a ripple of 0.1, which produces a filter of 764 taps.


> I got the result, but still not the same as the one from Matlab:
> 1. The channels seem to be shifted by 20 then wrapped around, i.e. assuming
> channel # starts with 1, the pfb result's channel 1 seems to be matlab's
> channel 21, and channel 20 in pfb result seems to be matlab's channel 40,
> channel 21 in pfb result seems to be matlab's channel 1.
> 2. Assuming the channel shift/wrap is done, the time-domain plot still are
> not the same, pfb channel 1 (matlab's channel 21)'s plot seems to be the
> same, and no signal channels do fit, but other channels' plot are not the
> same (i.e. the situation is similar to my result of DDC).

This strikes me as a simple naming issue; which channel you decide is
0 versus N. The way that I have the output of the channelizer working,
channel 0 is the DC channel, channel 1 is the next one to the "right"
(the next one in moving in the positive frequency direction). Channel
floor(N/2) will either span the last part of the positive spectrum and
the first part of the negative spectrum (if N is even) or be the last
channel in the positive spectrum (if N is odd). Channel N-1 is
directly "left" of channel 0.

It's hard to know what you mean by the signals being "the same." Are
they completely different?

> Is this because the filter taps in GNURadio is different from Matlab? Is
> there a way I can create a filter in GNURadio from a list of numbers
> exported from Matlab?

It could be a filter issue. If you want to, generate you prototype
filter in Matlab, then you can just copy-and-paste those into your GNU
Radio program as a list and pass that to the pfb_channelizer.

Tom