DSPRelated.com
Forums

Fir filter Hamming window 161 order bandpass implementation.

Started by vpot8787 6 years ago31 replieslatest reply 6 years ago235 views

I'm running the fir filter hamming window order 161 bandpass with pass frequency range 20hz-200hz,  but out put doesn't sound like the filter is working. it's just distorted output. following is the code and 162 coefficients. sampling rate 16khz. Can any body help why my filter not working?

coeff:-0.00035514,-0.0003342,-0.00031519,-0.0002975,-0.00028048,-0.00026346,-0.00024575,-0.00022664,-0.00020539,-0.00018126,-0.00015351,-0.00012138,-8.4119e-05,-4.0977e-05,8.7845e-06,6.5896e-05,0.00013107,0.00020499,0.00028834,0.00038174,0.00048581,0.00060111,0.00072817,0.00086749,0.0010195,0.0011846,0.0013631,0.0015553,0.0017614,0.0019817,0.0022161,0.0024647,0.0027276,0.0030046,0.0032955,0.0036002,0.0039183,0.0042495,0.0045934,0.0049494,0.005317,0.0056956,0.0060845,0.006483,0.0068902,0.0073053,0.0077274,0.0081555,0.0085887,0.0090258,0.0094659,0.0099077,0.01035,0.010792,0.011233,0.01167,0.012103,0.012531,0.012952,0.013366,0.013771,0.014165,0.014548,0.014918,0.015275,0.015617,0.015943,0.016253,0.016544,0.016817,0.017071,0.017304,0.017515,0.017705,0.017873,0.018018,0.018139,0.018237,0.01831,0.018359,0.018384,0.018384,0.018359,0.01831,0.018237,0.018139,0.018018,0.017873,0.017705,0.017515,0.017304,0.017071,0.016817,0.016544,0.016253,0.015943,0.015617,0.015275,0.014918,0.014548,0.014165,0.013771,0.013366,0.012952,0.012531,0.012103,0.01167,0.011233,0.010792,0.01035,0.0099077,0.0094659,0.0090258,0.0085887,0.0081555,0.0077274,0.0073053,0.0068902,0.006483,0.0060845,0.0056956,0.005317,0.0049494,0.0045934,0.0042495,0.0039183,0.0036002,0.0032955,0.0030046,0.0027276,0.0024647,0.0022161,0.0019817,0.0017614,0.0015553,0.0013631,0.0011846,0.0010195,0.00086749,0.00072817,0.00060111,0.00048581,0.00038174,0.00028834,0.00020499,0.00013107,6.5896e-05,8.7845e-06,-4.0977e-05,-8.4119e-05,-0.00012138,-0.00015351,-0.00018126,-0.00020539,-0.00022664,-0.00024575,-0.00026346,-0.00028048,-0.0002975,-0.00031519,-0.0003342,-0.00035514

code:-

for ( sec = 0 ; sec < 30 ; sec++ )
{
for ( msec = 0 ; msec < 1000 ; msec++ )
{
for ( sample = 0 ; sample < 16 ; sample++ )
{
/* Read 16-bit left channel Data */
I2S_readLeft(&data1);

/* Read 16-bit right channel Data */
I2S_readRight(&data2);


// As a new sample is received put it in two places
x[index1] = data1;
x[index2] = data1;
z[index1] = data2;
z[index2] = data2;
y1new=0;
y2new=0;
for(i=0;i<N;i++){
y1new = y1new + coeff[i]*x[index2-i];
y2new = y2new + coeff[i]*z[index2-i];
}

index1 = (index1+1)%N;
index2 = index1+N;
data1=y1new;
data2=y2new;
/* Write 16-bit left channel Data */
I2S_writeLeft(data1);

/* Write 16-bit right channel Data */
I2S_writeRight(data2);
}
if(sw3Pressed == TRUE)
{
break;
}

}
if(sw3Pressed == TRUE)
{
break;
}

}

[ - ]
Reply by BobLawlorApril 5, 2018

Hi,

I had a quick look at the frequency response of your filter using freqz() in Octave and it looks like a lowpass filter and not the bandpass one which you want.

You might have made a mistake in designing your filter. A good way to verify your filter design is to use freqz() to verify that the magnitude response looks right.

When it does I'd also recommend testing it with two pure synthetic tones (separately). One in the passband and one in the stopband.

Hope this helps

Bob


[ - ]
Reply by vpot8787April 5, 2018

Hi bob,

Thanks for your reply. I've used matlab fir1 function as shown in code below to find coefficients. I want to pass frequency between 20hz-200hz with 16khz sampling rate of codec.

fl=20*2/16000;

fh=200*2/16000;

wn=[fl,fh];

b=fir1(161,wn,'bandpass');

csvwrite('coeff.h',b);

[ - ]
Reply by BobLawlorApril 5, 2018

Hi,

I've just added this line of code to yours:

freqz(b,1,1024)

and it shows that the resulting filter isn't the BPF you're looking for.

Interestingly, if you push the filter order right up (I just tried 561 and 1061) you'll see that it is trying to give you the desired BPF.

I think the problem is that your desired Q-factor is so high that you need a much higher order filter.

It might be worth looking at using an IIR filter instead which can give high Q-factor and lower order. Just check the phase response in the passband as it might not be linear which can be a problem in some applications e.g. audio. 

Bob

[ - ]
Reply by vpot8787April 5, 2018

Hi bob,

561 is too high that my filter won’t be executed in real time. This is a digital stethoscope application and there shouldn’t be any noticeable latency

[ - ]
Reply by bmoersApril 5, 2018

The code does not show variable types or initialization values.  I assume your x and z buffers are sized accordingly (2N).  To my knowledge, the C5545 is a fixed point processor.  How are you handling the floating point coefficients given?  Presumably converted to 16 bit numbers.  Is the scaling of the multiplications handled correctly?

Others comments on the short filter length for the high degree of oversampling should also be considered in the design.


[ - ]
Reply by vpot8787April 5, 2018

bmoers,

thanks for the reply,

yes i didn't add the initialization part here.

well data is int16 data type, i'm just assigning the sum to data1 and data2. isn't that enough?

but if you plot the coeff with freqz it shows a good response, only it doesn't work in code 

[ - ]
Reply by bmoersApril 5, 2018

There are multiple problems with the implementation. Your need to scale the floating point coefficients so as to represent them as 16 bit integers.    Scaling so as to prevent overflow in fixed point is should not be trivialized.  Then y1new and y2new need to be 32 bit integers as the product of a 16 bits coefficient and a 16 bit data value results in 32 bits.  Also the sum of 162, 32 bit products must not overflow.  This is where the scaling of coefficients is important.  DSP hardware have wider accumulators to help prevent a sum overflow. The straight "C" implementation probably does not leverage such features.  Probably best to use a FIR library function crafted for your processor.  Finally you need to scale back your 32 bit sum of products to fit in the 16 bit output value.  Grab a book and start studying these concepts.

[ - ]
Reply by vpot8787April 5, 2018

here is the code for scaling of coeffs I did. But it still doesn't work

  int N=269;

    int i,n;

    int index1=0;

    int index2=N;

    float x[2*269];

    float z[2*269];

    float y1new,y2new;

    I2S_writeLeft(0);

I2S_writeRight(0);

for(n=0;n<N;n++){

coeff[n]=round(coeff[n]*256);

}

    for ( sec = 0 ; sec < 30 ; sec++ )

        {

            for ( msec = 0 ; msec < 1000 ; msec++ )

            {

                for ( sample = 0 ; sample < 16 ; sample++ )

                {

                    /* Read 16-bit left channel Data */

                I2S_readLeft(&data1);

                    /* Read 16-bit right channel Data */

                    I2S_readRight(&data2);

                    // As a new sample is received put it in two places

                    x[index1] = data1;

                    x[index2] = data1;

                    z[index1] = data2;

                    z[index2] = data2;

                    y1new=0;

                    y2new=0;

                    for(i=0;i<N;i++){

                    y1new = y1new + coeff[i]*x[index2-i];

                    y2new = y2new + coeff[i]*z[index2-i];

                    }

                    index1 = (index1+1)%N;

                    index2 = index1+N;

                            data1=round(y1new/256);

                            data2=round(y2new/256);

                    /* Write 16-bit left channel Data */

                    I2S_writeLeft(data1);

                    /* Write 16-bit right channel Data */

                    I2S_writeRight(data2);

                }

                if(sw3Pressed == TRUE)

                {

                break;

                }

            }

                if(sw3Pressed == TRUE)

                {

                break;

                }

         } 

[ - ]
Reply by dszaboApril 5, 2018

Oh snap!  Great catch.  If I’m not mistaken, all the coefficients are between -1 and 1, so the negative coefficients round to -1 and the positive coefficients round to zero.  This would almost certainly overflow the data during the summing operation, if not radically change the filter characteristics

Edit: It would also imply that the overall gain of the filter would be fairly large, have you observed that the volume of the audio is considerably louder with the filter running than without

[ - ]
Reply by vpot8787April 5, 2018

Volume decreases and distorts when filter is run,

I've used the above scaling technique for coeffs. It still distorts, doesn't sound like the filter should actually sound.

[ - ]
Reply by dszaboApril 5, 2018

Your code is very clearly using floating point math, so scaling the coefficients won’t change the results in a meaningful way.  You should really measure the execution time for the filter.  If that device does not have a FPU then in all likelihood, this code is taking waaaaaay longer to execute than you thought.

[ - ]
Reply by vpot8787April 5, 2018

dszabo,

i'm using TMS320C5545 boosterpack dsp which is Fixed-Point Digital Signal Processor. I'm pretty new to dsp n digital filters, but I tried to convert the floating coeff to 16bit integer as show in my updated code below. can you point me where i'm going wrong.sampling rate of codec set to 8000hz

Int16 coeff1[NUMCOEFFS];

Int16 data1,data2;

     int N=17;

    int i,n;

    int index1=0;

    int index2=N;

    Int16 x[2*17];

    Int16 z[2*17];

    float y1new,y2new;

    I2S_writeLeft(0);

I2S_writeRight(0);

for(n=0;n;n++){>

coeff1[n]=round(coeff[n]*65536);

}

    for ( sec = 0 ; sec < 30 ; sec++ )

        {

            for ( msec = 0 ; msec < 1000 ; msec++ )

            {

                for ( sample = 0 ; sample < 8 ; sample++ )

                {

                    /* Read 16-bit left channel Data */

                I2S_readLeft(&data1);

                    /* Read 16-bit right channel Data */

                    I2S_readRight(&data2);

                    // As a new sample is received put it in two places

                    x[index1] = data1;

                    x[index2] = data1;

                    z[index1] = data2;

                    z[index2] = data2;

                    y1new=0;

                    y2new=0;

                    for(i=0;i;i++){>

                    y1new = y1new + coeff1[i]*x[index2-i];

                    y2new = y2new + coeff1[i]*z[index2-i];

                    }

                    index1 = (index1+1)%N;

                    index2 = index1+N;

                            data1=round(y1new/65536);

                            data2=round(y2new/65536);

                    /* Write 16-bit left channel Data */

                    I2S_writeLeft(data1);

                    /* Write 16-bit right channel Data */

                    I2S_writeRight(data2);

                }

                if(sw3Pressed == TRUE)

                {

                break;

                }

            }

                if(sw3Pressed == TRUE)

                {

                break;

                }

        }

[ - ]
Reply by dszaboApril 5, 2018

A couple things. First your accumulators are still floating point, change them to 32 bit integers. Second, you are using 16 bits of precision for your coefficients.  As already pointed out, you should be cautious of overflow, either by reduction the precision of the coefficients or leveraging the instruction set of the device. Third, it looks like a couple of your for loops are declared incorrectly

[ - ]
Reply by vpot8787April 5, 2018

dszabo,

Thanks for the suggestions. I have made some great progress with the following changes in code. I main issue was over flow so then I reduced precision from 256 to 32. I've recorded a despacito song filter with passband 20-200hz with hamming window 161 order FIR. output is lil distorting at high beats. attached is audio link ,listen to it. There is a hum in the recording that's from laptop adc, output from dsp doesn't has hum, only distortion. Also attached spectrum plot of recording.

Recorded_Audio (Listen audio)

Plot

[ - ]
Reply by dszaboApril 5, 2018

It sounds to me like you are overflowing your output. Modern music is mastered to be really loud, and if you look closely at some waveforms, you’ll find instances of what looks kind of like clipping. Clipped waveforms have more high frequency content then non clipped. If I filter out some of the frequency content that was created by clipping, the waveform will extend beyond the clipping point. 

The point is: just because your input signal is 16-bit and your filter does not exceed unity gain, that does not generally mean your output is 16bit as well, even though it probably is most of the time. 

If it were me, I’d probably add a saturation stage right before the output to limit values to 16-bit, because I’d rather have saturation distortion than overflow.  In any case, you’ll want to turn down the input a bit.

You can probably do better on your coefficient precision then, at least put it back up to 8 bit. You can also reduce the precision of your input waveform.  You’ll end making some kind of trade based on what you think is best

[ - ]
Reply by vpot8787April 5, 2018

dszabo,

1)For some reason filter is passing frequencies upto 4000hz or so, I've designed the filter coeff for bandpass of 20-200hz. Please look at the coeff i used below. if you plot it look perfect banpass in 20-200hz, but when implemented in code it doesn't sound like it should.you can hear the clip in my previous reply. i also imported the original sound track into audacity and did a lowpass with cutoff 200 and roll-off 24db. That sounded right(it should sound like if you turn on loud music with good bass in your car and close the doors and listen it from outside near to car (for a passband of 20-200hz). but my filter implementation sounds like close to original audio;meaning able to hear high frequencies, but with low volume kindof affect allover the time. 

2)I was speaking to a TI specialist in general about the high order(161 order) filter implementation and distortion issues I faced before. He suggested me to use Q15 format. Both input and coeffs data being implemented should be presented in Q15 format.see below for his reply.

a) Yes. C5545 is a fixed point DSP, but the way we do the filtering is using the Q15 (1 bit sign and 15 bit data, the value is between -1 and 1, i.e. the decimal point is after bit 15). Both data and coefficients are presented in Q15 format. To do the filtering algorithm, you can do the following:long res_q30, res_q15;int data_q15, coef_q15;res_q30 = (long)data_q15*coef_q15;res_q15 = res_q30>>15;

b) Most likely is the incorrect filtering algorithm or the not using Ping Pong buffers (one is used for receiving data, the other is used for data processing). You can use the Audio codec loopback example as the start point (ccs_v6.x_examples\i2s\CSL_I2S_AudioCodec_DMA), then add your filtering algorithm to UserAlgorithm() in AudioCodec_DMA.c.You can also download the following audio filtering demo for C5515 eZDSP Kit from https://code.google.com/archive/p/c5505-ezdsp/downloadshttps://storage.googleapis.com/google-code-archive-downloads/v2/code.google.com/c5505-ezdsp/C5515_eZdsp_Audio_Filter_Demo1.zipThe filtering algorithm is implemented in assembly (fir() in fir.asm). The main loop is in main_bypass1.c. The I2S used in C5515 ezdsp is different than the one used in C5545 ezdsp though.Both examples are using Ping-Pong Buffers.

coeff:

0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.0001,-0.0001,-0.0001,-0.0001,-0.0001,-0.0002,-0.0002,-0.0001,-0.0001,-0.0001,-0.0001,0,0,0.0001,0.0002,0.0002,0.0003,0.0004,0.0004,0.0004,0.0004,0.0004,0.0004,0.0003,0.0002,0.0001,0,-0.0002,-0.0003,-0.0005,-0.0006,-0.0008,-0.0009,-0.0009,-0.001,-0.001,-0.0009,-0.0008,-0.0006,-0.0003,-0.0001,0.0002,0.0005,0.0009,0.0012,0.0015,0.0017,0.0019,0.0019,0.0019,0.0018,0.0016,0.0013,0.0009,0.0004,-0.0001,-0.0007,-0.0013,-0.0019,-0.0025,-0.003,-0.0034,-0.0036,-0.0037,-0.0036,-0.0033,-0.0028,-0.0022,-0.0013,-0.0004,0.0007,0.0018,0.0029,0.004,0.005,0.0058,0.0064,0.0067,0.0067,0.0064,0.0058,0.0047,0.0034,0.0017,-0.0002,-0.0023,-0.0044,-0.0066,-0.0087,-0.0105,-0.012,-0.0131,-0.0136,-0.0135,-0.0128,-0.0112,-0.0089,-0.0058,-0.002,0.0026,0.0077,0.0135,0.0196,0.0259,0.0324,0.0388,0.0449,0.0507,0.0559,0.0603,0.0639,0.0666,0.0682,0.0687,0.0682,0.0666,0.0639,0.0603,0.0559,0.0507,0.0449,0.0388,0.0324,0.0259,0.0196,0.0135,0.0077,0.0026,-0.002,-0.0058,-0.0089,-0.0112,-0.0128,-0.0135,-0.0136,-0.0131,-0.012,-0.0105,-0.0087,-0.0066,-0.0044,-0.0023,-0.0002,0.0017,0.0034,0.0047,0.0058,0.0064,0.0067,0.0067,0.0064,0.0058,0.005,0.004,0.0029,0.0018,0.0007,-0.0004,-0.0013,-0.0022,-0.0028,-0.0033,-0.0036,-0.0037,-0.0036,-0.0034,-0.003,-0.0025,-0.0019,-0.0013,-0.0007,-0.0001,0.0004,0.0009,0.0013,0.0016,0.0018,0.0019,0.0019,0.0019,0.0017,0.0015,0.0012,0.0009,0.0005,0.0002,-0.0001,-0.0003,-0.0006,-0.0008,-0.0009,-0.001,-0.001,-0.0009,-0.0009,-0.0008,-0.0006,-0.0005,-0.0003,-0.0002,0,0.0001,0.0002,0.0003,0.0004,0.0004,0.0004,0.0004,0.0004,0.0004,0.0003,0.0002,0.0002,0.0001,0,0,-0.0001,-0.0001,-0.0001,-0.0001,-0.0002,-0.0002,-0.0001,-0.0001,-0.0001,-0.0001,-0.0001,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
[ - ]
Reply by vpot8787April 5, 2018

Update:-

filter is definitely doing something, but for some reason its folding high frequencies. I've designed a lowpass filter with fc=200 and fstop=500 and sampling rate=8000hz, and order 161. see attached matlab pic.

filter_design_9287.gif

using those coeffs I ran filter and following is the frequency plot. see how frequencies after 500 are passing and also getting fold.

dsp_filtered_58534.gif

and below is the frequency plot for original audio song 

original_21405.gif

below is after applying audacity lowpass filter for cutoff 200hz and roll off 24db, for above original audio song. I would like to achieve such output with my dsp filter.

audacity_filtered(24db)_67883.gif

[ - ]
Reply by dszaboApril 5, 2018

I'm sorry to hear that you are still having troubles with your implementation.

If you want to get your code working, I would strongly recommend creating a version of your code that runs on your PC such that your read data from an audio file and output to a .WAV file.  Use that to create a filter that works.

In terms of your target application, just listen to your specialist and plug your desired coefficients into his/her examples.  They are going to be much better optimized than anything you'll ever write in C.

Best of luck.

[ - ]
Reply by vpot8787April 5, 2018

hi dszabo,

If you want to get your code working, I would strongly recommend creating a version of your code that runs on your PC such that your read data from an audio file and output to a .WAV file.  Use that to create a filter that works.

do you mean to say run a code on dsp and ccs to read the data from despacito audio file i was using and output that to a .wav file?and the use the .wav file as input for filtering code? what difference would it make?

Thanks.


[ - ]
Reply by dszaboApril 5, 2018

That is not what I meant to say. Sorry I wasn’t able to help more effectively

[ - ]
Reply by vpot8787April 5, 2018

I really appreciate your advise so far.

what code were you talking about ,reading data from audio file and output to a .wav file?

[ - ]
Reply by dszaboApril 5, 2018

I’ll try my best to phrase it a little differently. Create a command line tool or other application on your PC that opens a wav file, processes it with your filter, then puts the result in an output wav file. This removes filter from both the target device and the complexities of real time processing. You can then use audacity as you already have to verify that result matches your expectation. This will allow you to determine whether there is an issue with your filter implementation or with how you are integrating it into the target.  You can also debug the code without worrying about real time data and remove any alterations resulting from the ADC and DAC processes.

[ - ]
Reply by dszaboApril 5, 2018

I deleted my last post because I didn’t understand that you were duplicating you delay buffer to get around the modulo addressing problem. That’s really clever.

I assume that when you are saying that the output of the filter is distorted, you are referring to your target system running in real time. I can’t find any fault with your filter implementation, aside from the fact that you may not be using enough coefficients to get the filter performance you require, but that shouldn’t cause any distortion.

Have you ruled out that the distortion may be a result of your I2S peripheral?  I2S peripherals can be a little finicky if you aren’t dealing with the buffering properly. What is your target MCU?  Also, have you made sure that your filter is executing fast enough to keep up with your sample rate?

If you want to validate your filter implementation, I would recommend breaking the code into its own module and testing it outside the target.

[ - ]
Reply by vpot8787April 5, 2018

Hi dszabo

I’m using c5545 booster pack dsp which comes inbuilt with aic3206

I’ve set codec sampling rate as 16khz

I’ve tested audio line in loop back test and it outputs aydio without distorting. I thought distortion was due to the filter coefficients. I passed a signal tone of my pass band frequency and it doesn’t work or it even passes frequencies even more than 5khz. I’m interested in bandpass of 20-200h

[ - ]
Reply by dszaboApril 5, 2018

It may be that the latency introduced by the filter execution is causing your transmit buffer to run dry.  On the cc5545, you could check the OUERRFL bit in the I2SINTFL register to make sure no over/under runs occurred.

You’ll also want to be sure that you aren’t saturating the output.  Does the distortion happen with low level signals or just high level signals

[ - ]
Reply by vpot8787April 5, 2018

Hi dszabo,

I’m sure it’s not output saturation because audio loopback tst worked without distortion. Anyways I’m playing a audio song file as input from laptop

Can you specify what is the OUERRFL bit value I should set. Thanks in ad

[ - ]
Reply by dszaboApril 5, 2018

It’s an interrupt status bit, you can read it to check for errors. Refer to the MCU user guide for more detail.

You might also try transmitting a zero before you go into your loop, this might keep the buffer from running dry, assuming that is what the problem is

[ - ]
Reply by vpot8787April 5, 2018

dszabo,

you mean before the following code you want me to give, I2S_writeLeft(0) and I2S_writeRight(0)??

code:-

for ( sec = 0 ; sec < 30 ; sec++ )
{
for ( msec = 0 ; msec < 1000 ; msec++ )
{
for ( sample = 0 ; sample < 16 ; sample++ )
{
/* Read 16-bit left channel Data */
I2S_readLeft(&data1);

/* Read 16-bit right channel Data */
I2S_readRight(&data2);


[ - ]
Reply by dszaboApril 5, 2018

Pretty much yes, it’s a pretty easy test without having to overly complicate your implementation. Ideally, you’d be managing a double buffer and using DMA, but that’s a big change if you aren’t doing it already

[ - ]
Reply by vpot8787April 5, 2018

i tried it but it didn't work :/