jgaeddert / liquid-dsp

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

Assert when designing bandpass #122

Open viraptor opened 6 years ago

viraptor commented 6 years ago

I've tried to use a bandpass based on the example in http://liquidsdr.org/doc/firdes/

Unfortunately, I get result:

src/filter/src/firdespm.c:711: firdespm_iext_search: Assertion `num_found < nmax' failed.

I'm not sure what it means exactly and what's the way to deal with this? I can't find anything related in the docs. I'm using the LIQUID_FIRDESPM_BANDPASS type and the definition:

    float bands[6] = {
        0.0f, 18.5e3/sample_rate,
        18.7e3/sample_rate, 19.3e3/sample_rate, 
        19.5e3/sample_rate, 0.5f,
    };
    float des[3] = { 0.0f, 1.0f, 0.0f };
    float weights[3] = { 1.0f, 1.0f, 1.0f };

    liquid_firdespm_wtype wtype[3] = {
        LIQUID_FIRDESPM_EXPWEIGHT,
        LIQUID_FIRDESPM_EXPWEIGHT,
        LIQUID_FIRDESPM_EXPWEIGHT,
    };
jgaeddert commented 6 years ago

Unfortunately this is an internal error. The firdespm method is an implementation of the Park-McClellan algorithm for minimax FIR filter design that iterates over the Remez exchange algorithm. Theoretically there are supposed to be a certain number of extremal frequencies, but sometimes (because of machine precision) more are found.

Right now I haver an assertion that the number is valid, but this will kick the program out rather than exiting gracefully. I should try to be more resilient and/or perhaps a bit nicer if it encounters such an error.

viraptor commented 6 years ago

To be more specific, I didn't actually mind the assert (but sure - an error would be nicer). What I was confused about was - so what do I do about it?

Would it be possible to add some suggestions? I removed the wtype and replaced it with null like I saw in some example and it suddenly started working. But I'm left with "huh?" :-)

jgaeddert commented 6 years ago

Well, there are a few things I can suggest for the time being. To begin, you define three regions without much guard in between them (originally I think you didn't have any guard band in between them); increasing that might help. The algorithm is pretty sensitive to the parameters of the filter design, but I might be able to improve it if I could reproduce your results. I've guessed at your parameters (sample_rate = 100e3f;) and I can verify on my end that I'm getting the same error.

jgaeddert commented 6 years ago

As an aside, if you just want to create a simple band-pass filter without too much complexity, you can do so with a low-pass windowed-Kaiser and then frequency shifting it. I've attached an example to do just that. I know this isn't exactly what your original ticket was about (firdespm), but it might solve your larger problem. firdeswk

#include <stdio.h>
#include <math.h>
#include <liquid/liquid.h>

int main() {
    // define filter length, type, number of bands
    unsigned int n           =  251;        // filter length
    float        freq_center =  20e3f;      // center frequency [Hz]
    float        bandwidth   =   2.0e3f;    // filter bandwidth [Hz]
    float        sample_rate = 100e3f;      // sample rate [Hz]
    float        As          =  40.0f;      // stop-band suppression [dB]

    // allocate memory for array and design low-pass prototype filter
    float h[n];
    liquid_firdes_kaiser(n, 0.5f*bandwidth/sample_rate, As, 0.0f, h);

    // shift frequency
    unsigned int i;
    for (i=0; i<n; i++)
        h[i] *= cosf(2*M_PI*((float)i - 0.5f*(float)n) * freq_center / sample_rate);

    // generate output file for plotting
    FILE * fid = fopen("firdeswk.m","w");
    fprintf(fid,"clear all; close all;\n");
    fprintf(fid,"nfft = 1200;\n");
    for (i=0; i<n; i++)
        fprintf(fid,"h(%3u) = %12.8f;\n", i+1, h[i]);
    fprintf(fid,"f = ([0:(nfft-1)]/nfft-0.5)*%12.4e;\n", sample_rate*1e-3f);
    fprintf(fid,"plot(f,20*log10(abs(fftshift(fft(h,nfft)))));\n");
    fprintf(fid,"grid on;\n");
    fprintf(fid,"xlabel('Frequency [kHz]');\n");
    fprintf(fid,"ylabel('Power Spectral Density [dB]');\n");
    fclose(fid);
}
viraptor commented 6 years ago

You can shift low pass to band pass by multiplying taps? Whaaaaaaat? I learned something cool today, thanks!

fsheikh commented 6 years ago

On an intel core i7 machine I was not able to reproduce the problem. And octave was also able to design the filter with 251 filter length. Frequency response from liquid-dsp with these parameters. firdespm_bandpass_example