1
$\begingroup$

I want to filter out complex valued samples - Bandpass filter. I have 100 complex samples and I wish to filter out a range of frequencies from it, as you can see in the result of the FFT, there are 3 frequencies here.

  1. First At 252KHz
  2. Second At 315KHz
  3. Third At 378KHz

enter image description here I want to only keep the middle one (315KHz)

I am using this simple FIR filter class by by Mike Perkins class, it already includes 3 filters, LPF, HPF and BPF which is what I need.

The problem is, it was designed to work with real samples, i.e wav files, now I am trying to get around this but I'm having troubles doing so, you can find documentation of his class in the github link, but I'm going to post the filter class here below,.

I'm going to focus on the do_sample(float data_sample) function, this function takes a real data sample as input, and it returns real filtered data sample.

Here is the code for it:

float Filter::do_sample(float data_sample)
{
    int i;
    float result;

    if( m_error_flag != 0 ) return(0);

    for(i = m_num_taps - 1; i >= 1; i--){
        m_sr[i] = m_sr[i-1];
    }   
    m_sr[0] = data_sample;

    result = 0;
    for(i = 0; i < m_num_taps; i++) result += m_sr[i] * m_taps[i];

    return result;
}

What I've tried to do:

  1. I tried sending the real and imaginary parts separately to the function
  2. I tried making 2 separate Filter objects, one for the real part and the other for the imaginary parts and using the function with each object

both "solutions" didn't work and returning me the same result:

enter image description here

Finally, here is the filter class:

filt.h

#ifndef _FILTER_H
#define _FILTER_H

#define MAX_NUM_FILTER_TAPS 1000
#pragma warning(disable:4996)

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <inttypes.h>
#include <iostream>

enum filterType {LPF, HPF, BPF};

class Filter{
    private:
        filterType m_filt_t;
        int m_num_taps;
        int m_error_flag;
        float m_Fs;
        float m_Fx;
        float m_lambda;
        float *m_taps;
        float *m_sr;
        void designLPF();
        void designHPF();

        // Only needed for the bandpass filter case
        float m_Fu, m_phi;
        void designBPF();

    public:
        Filter(filterType filt_t, int num_taps, float Fs, float Fx);
        Filter(filterType filt_t, int num_taps, float Fs, float Fl, float Fu);
        ~Filter();
        void init();
        float do_sample(float data_sample);
        int get_error_flag(){return m_error_flag;};
        void get_taps( float *taps );
};

#endif

filt.cpp:

#include "filt.h"

#define ECODE(x) {m_error_flag = x; return;}

// Handles LPF and HPF case
Filter::Filter(filterType filt_t, int num_taps, float Fs, float Fx)
{
    m_error_flag = 0;
    m_filt_t = filt_t;
    m_num_taps = num_taps;
    m_Fs = Fs;
    m_Fx = Fx;
    m_lambda = M_PI * Fx / (Fs/2);

    if( Fs <= 0 ) ECODE(-1);
    if( Fx <= 0 || Fx >= Fs/2 ) ECODE(-2);
    if( m_num_taps <= 0 || m_num_taps > MAX_NUM_FILTER_TAPS ) ECODE(-3);

    m_taps = m_sr = NULL;
    m_taps = (float*)malloc( m_num_taps * sizeof(float) );
    m_sr = (float*)malloc( m_num_taps * sizeof(float) );
    if( m_taps == NULL || m_sr == NULL ) ECODE(-4);
    
    init();

    if( m_filt_t == LPF ) designLPF();
    else if( m_filt_t == HPF ) designHPF();
    else ECODE(-5);

    return;
}

// Handles BPF case
Filter::Filter(filterType filt_t, int num_taps, float Fs, float Fl,
               float Fu)
{
    m_error_flag = 0;
    m_filt_t = filt_t;
    m_num_taps = num_taps;
    m_Fs = Fs;
    m_Fx = Fl;
    m_Fu = Fu;
    m_lambda = M_PI * Fl / (Fs/2);
    m_phi = M_PI * Fu / (Fs/2);

    if( Fs <= 0 ) ECODE(-10);
    if( Fl >= Fu ) ECODE(-11);
    if( Fl <= 0 || Fl >= Fs/2 ) ECODE(-12);
    if( Fu <= 0 || Fu >= Fs/2 ) ECODE(-13);
    if( m_num_taps <= 0 || m_num_taps > MAX_NUM_FILTER_TAPS ) ECODE(-14);

    m_taps = m_sr = NULL;
    m_taps = (float*)malloc( m_num_taps * sizeof(float) );
    m_sr = (float*)malloc( m_num_taps * sizeof(float) );
    if( m_taps == NULL || m_sr == NULL ) ECODE(-15);
    
    init();

    if( m_filt_t == BPF ) designBPF();
    else ECODE(-16);

    return;
}

Filter::~Filter()
{
    if( m_taps != NULL ) free( m_taps );
    if( m_sr != NULL ) free( m_sr );
}

void 
Filter::designLPF()
{
    int n;
    float mm;

    for(n = 0; n < m_num_taps; n++){
        mm = n - (m_num_taps - 1.0) / 2.0;
        if( mm == 0.0 ) m_taps[n] = m_lambda / M_PI;
        else m_taps[n] = sin( mm * m_lambda ) / (mm * M_PI);
    }

    return;
}

void 
Filter::designHPF()
{
    int n;
    float mm;

    for(n = 0; n < m_num_taps; n++){
        mm = n - (m_num_taps - 1.0) / 2.0;
        if( mm == 0.0 ) m_taps[n] = 1.0 - m_lambda / M_PI;
        else m_taps[n] = -sin( mm * m_lambda ) / (mm * M_PI);
    }

    return;
}

void 
Filter::designBPF()
{
    int n;
    float mm;

    for(n = 0; n < m_num_taps; n++){
        mm = n - (m_num_taps - 1.0) / 2.0;
        if( mm == 0.0 ) m_taps[n] = (m_phi - m_lambda) / M_PI;
        else m_taps[n] = (   sin( mm * m_phi ) -
                             sin( mm * m_lambda )   ) / (mm * M_PI);
    }

    return;
}

void 
Filter::get_taps(float* taps)
{
    int i;

    if (m_error_flag != 0) return;

    for (i = 0; i < m_num_taps; i++) taps[i] = m_taps[i];

    return;
}

void 
Filter::init()
{
    int i;

    if( m_error_flag != 0 ) return;

    for(i = 0; i < m_num_taps; i++) m_sr[i] = 0;

    return;
}

float 
Filter::do_sample(float data_sample)
{
    int i;
    float result;

    if( m_error_flag != 0 ) return(0);

    for(i = m_num_taps - 1; i >= 1; i--){
        m_sr[i] = m_sr[i-1];
    }   
    m_sr[0] = data_sample;

    result = 0;
    for(i = 0; i < m_num_taps; i++) result += m_sr[i] * m_taps[i];

    return result;
}

UPDATE: I added a main and recreated this issue so you can run it for yourself.

main.cpp:

#include <vector>
#include "filt.h"
#include <iostream>
#include <complex>

using std::vector;
using std::complex;

typedef vector<complex<double>> complexSignal;

vector<double> doFilter(Filter* my_filter, const complexSignal& samples, int N, bool real_part)
{
    vector<double> to_return;
    float part;
    for (int i = 0; i < N; i++)
    {
        if (real_part)
            part = my_filter->do_sample(samples[i].real());
        else
            part = my_filter->do_sample(samples[i].imag());
        to_return.push_back(part);
    }
    return to_return;

}

int main(int argc, char *argv[])
{
    complexSignal complex_samples;
    int frequency1 = 252000;
    int frequency2 = 315000;
    int frequency3 = 378000;

    double Fs = 1260000.0;
    const double  T = 1 / Fs;

    //generate complex sine wave
    for (int i = 0; i < 100; i++)
    {
        complex_samples.push_back({ (0.7 * cos(2 * M_PI * frequency1 * (i * T))), (0.7 * sin(2 * M_PI * frequency1 * (i * T))) });
        complex_samples[i] += { (0.7 * cos(2 * M_PI * frequency2 * (i * T))), (0.7 * sin(2 * M_PI * frequency2 * (i * T))) };
        complex_samples[i] += { (0.7 * cos(2 * M_PI * frequency3 * (i * T))), (0.7 * sin(2 * M_PI * frequency3 * (i * T))) };
    }

    Filter* BPF_Filter_real;
    Filter* BPF_Filter_imag;

    BPF_Filter_real = new Filter(BPF, 51, 1260.0, 308.7, 321.3);
    BPF_Filter_imag = new Filter(BPF, 51, 1260.0, 308.7, 321.3);

    vector<double> real_filtered = doFilter(BPF_Filter_real, complex_samples, 100, 1);
    vector<double> imag_filtered = doFilter(BPF_Filter_imag, complex_samples, 100, 0);

    complexSignal complex_filtered;
    for (int i = 0; i < 100; i++)
    {
        complex_filtered.push_back({ real_filtered[i], imag_filtered[i] });
    }
}   

UPDATE 2: RESULTS: enter image description here

$\endgroup$
14
  • 1
    $\begingroup$ Filtering is linear, so $\textrm{BPF}\lbrace a(t)+jb(t) \rbrace = \textrm{BPF}\lbrace a(t) \rbrace + j \textrm{BPF}\lbrace b(t) \rbrace$. $\endgroup$
    – MBaz
    Dec 3 '21 at 17:22
  • $\begingroup$ So you're saying I should filter the real and imag part separately, but I already tried that (look in the post) @MBaz $\endgroup$ Dec 3 '21 at 17:47
  • 1
    $\begingroup$ Yeah -- but if you know it should work, the conclusion is that the problem is in your code, not in your approach. $\endgroup$
    – MBaz
    Dec 3 '21 at 18:24
  • $\begingroup$ I understand, so I created a small main file with how I implemented this and added it to the post, kindly take a look at let me know if you spot anything wrong, you can also run it for yourself. @MBaz $\endgroup$ Dec 3 '21 at 18:55
  • $\begingroup$ @yarinCohen you're simply not doing BPF(real) and BPF(imag); you're doing the identical call for both BPF_Filter_real and _imag. This is really just a debugging, not a signal processing problem, and using your debugger would have led you to the same conclusion as I just did about you calling doFilter with exactly the same arguments. $\endgroup$ Dec 3 '21 at 19:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.