jgaeddert / liquid-dsp

digital signal processing library for software-defined radios
http://liquidsdr.org
MIT License
1.89k stars 443 forks source link

Coefficients of msresamp_crcf #294

Closed guilleortas closed 2 years ago

guilleortas commented 2 years ago

Hi there, I'd like to know what are the coefficients used by the multi-stage resampler. I couldn't find information relative to the coefficients in the documentation. Thank you for this amazing library! Regards

jgaeddert commented 2 years ago

So there are multiple stages of resampling (hence its name). For interpolation, it's an arbitrary-rate resampler followed by a series of half-band resamplers. At current, there isn't an interface for getting the coefficients for each stage. What are you looking for, specifically?

guilleortas commented 2 years ago

Thank you for your reply. I was looking for a way of determining the filter used, as I'm getting very weird outputs with the 'resamp' and 'msresamp' objects. To illustrate, I'm attaching the spectrum of the input signal and the outputs I get with 'msresmp' as well as with the default parameters (using 'resamp_crcf_create_default') with 'resamp', respectively. In addition, for reference, I'm attaching the output of the rational resampler block of GNURadio, showing the expected result. In all cases I'm resampling the same 2ms segment of an input signal from 60MHz to 80MHz.

liquid-dsp_resamp_spectra gnuradio_fractional_resampler_spectrum

jgaeddert commented 2 years ago

Interesting. Can you provide the code snippet for how you're creating the object and how you're using it?

Alternatively you can use the rresamp_crcf object which implements a rational rate resampler. As per your example:

unsigned int P = 80; // interpolation factor
unsigned int Q = 60; // decimation factor
unsigned int m = 12; // filter semi-length
float bw = 0.5f; // resampling filter bandwidth
float As = 60.0f; // resampling filter stop-band suppression [dB]

rresamp_crcf q = rresamp_crcf_create_kaiser(P,Q,m,bw,As);    

float complex buf_0[Q]; // input
float complex buf_1[P]; // output

{
    // repeat as necessary
    rresamp_crcf_execute(q, buf_0, buf_1);
}

rresamp_crcf_destroy(q);
jgaeddert commented 2 years ago

Useful example here: https://github.com/jgaeddert/liquid-dsp/blob/master/examples/rresamp_crcf_example.c

guilleortas commented 2 years ago

Of course, here it is:

float resamp_ratio = rate_out/rate_in; // rate_in and rate_out defined as user inputs
float filt_bw = std::min((float)0.49, (float)0.5*resamp_ratio);
float As = 60;
unsigned int filt_semi_len = 4;
unsigned int filt_bank_size = 8;
resampler = resamp_crcf_create(resamp_ratio, filt_semi_len, filt_bw, As, filt_bank_size);
...
size_t samps_per_buffer = 128000;
std::vector<std::complex<float>> samples_in(samps_per_buffer);
unsigned int samples_out_size = ceilf(samps_per_buffer * resamp_ratio);
std::vector<std::complex<float>> samples_out(samples_out_size);
// samples_in is populated from an IQ file in every iteration
// samples_out is dumped to a file in every iteration, writing exactly 'samples_out_size' as returned by the execution of the resamp object 
{
    resamp_crcf_execute_block(resampler , &samples_in.front(), samps_per_buffer, &samples_out.front(), &samples_out_size);
}

I'll try the rresamp_crcf object and report back. Thank you again!

guilleortas commented 2 years ago

In the meantime, I can leave here a few output spectrums I got while trying to figure out what's going on when tweaking the resampler parameters. It's worth noting the shocking difference between setting a filter bandwidth of 0.499999 vs 0.49 leaving everything else the same. filter_out_comparison

jgaeddert commented 2 years ago

Right. You need to provision at least some amount of guard band otherwise the filter design becomes extremely difficult. I think that might be part of the issue here. The msresamp is trying to compute the minimum filter length for each stage given the constraints. With no guard on the edge of the band, that's a difficult problem.

I agree though, that the response is odd, even still.

guilleortas commented 2 years ago

I completely agree, these values are unnecessarily close to 0.5. I just used them for some testing I did for myself trying to figure out what was happening. Still, setting the filter bandwidth to 0.45 or 0.40 still yielded very weird results with mangled frequency components. I was looking at rresamp_crcf and I'm wondering if there's any way of running the resampler in lager batches than just Q input samples at a time. Can I somehow feed a a larger signal buffer like for msresamp or resamp?

jgaeddert commented 2 years ago

For resamp_crcf, you should probably increase the filter length and number of filters in the bank. Maybe use filter_semi_len = 12 and filt_bank_size = 64.

Yes, you can run rresamp_crcf in blocks. Amending my previous example, I think this should work fine:

unsigned int P = 80; // interpolation factor
unsigned int Q = 60; // decimation factor
unsigned int m = 12; // filter semi-length
float bw = 0.5f; // resampling filter bandwidth
float As = 60.0f; // resampling filter stop-band suppression [dB]

rresamp_crcf q = rresamp_crcf_create_kaiser(P,Q,m,bw,As);    

unsigned int num_blocks = 100;
float complex buf_0[100*Q]; // input
float complex buf_1[100*P]; // output

{
    // repeat as necessary
    //rresamp_crcf_execute(q, buf_0, buf_1);
    rresamp_crcf_execute_block(q, buf_0, num_blocks, buf_1);
}

rresamp_crcf_destroy(q);

Note that the true resampling rate is 80/60 = 4/3. Internally, the rational resampler uses this prime ratio but respects the original values of P and Q when you run it.

jgaeddert commented 2 years ago

See also: https://liquidsdr.org/api/rresamp_crcf/#execute_block

(Note that the API definition pages are a work in progress)

guilleortas commented 2 years ago

Hello again, thank you for your reply. The execution of the resampler in blocks is great. Regarding the original problem, after implementing my resampler using rresamp_crcf and the problem persisting, I was able to dig deeper to come to the realization that this behavior was purely caused by a mistake from my side when converting the result of the resampling to the desired output format. It was a corner case involving a type cast causing trouble with specific input signals which I haven't seen before, totally unconnected from liquid-dsp. I'm closing the issue, I was able to confirm that the resamplers msresamp_crcf, resamp_crcf and rresamp_crcf are working great, as expected. Thank you again for you attention.

jgaeddert commented 2 years ago

Glad to hear it worked out! Good luck!