3
$\begingroup$

I used PyFDA to design a low-pass FIR filter for the final stage of decimation of a complex (I/Q) stream, from 2048 MHz to 16kHz. The last stage, 32kHz to 16kHz is the only non half-band filter, as the transition needs to be wholly within the lower half of the frequency range to prevent aliasing. I'm using 14-bit data so I went for 90dB stopband attenuation. An FFT on the coefficients mostly concurs with the output of PyFDA. I am not scaling the FFT output at all.

LP FIR Order 214

To test the filter, I created a series of half-second complex sinusoids from -14.75 to +15.25 kHz stepping 1 kHz at a time with a peak amplitude of 0.9. Here's the signal's FFT before filtering. enter image description here

Here it is after filtering but before decimation. Filtering was done using np.convolve, but I get the same result with the production C code I'm working on. All is single precision floating point.

enter image description here

Just what I'd expect, although the intricacies of FFT scaling are still eluding me, which I believe might have something to do the problem.

Now the part that isn't working. I take the signal and convolve it with the filter. If I plot the RMS of the filtered waveform I don't see the attenuation that the filter is supposed to deliver. Rather than 90dB stopband attenuation I'm getting 60dB. Here are the RMS amplitudes for each 1/2 second interval.

enter image description here

I removed the beginning and end of each segment because the spikes were making it look even worse (about 10dB worse). This is the code that plotted RMS called with dwell=0.

    def plot_rms(data, sr, dwell=0):
        if 0 == dwell:
            dwell = sr / 2
    
        chunks = data.reshape(-1, dwell)[:, dwell // 16: dwell * 15 // 16]  
        rms = np.sqrt(np.mean(np.abs(chunks * chunks), axis=1))
        x_axis = np.linspace(0, data.size / sr, rms.size)
        plt.plot(x_axis, 20 * np.log10(rms))
        plt.show()

The filtered RMS also looks a bit weird. Zooming in on the stopband, I see this in the first 1.5 seconds. I can't imagine why it increases in intervals like that.

![enter image description here

So I'm definitely confused. How do I compare the designed & calculated frequency response with the actual RMS attenuation? Do I need to design a bigger filter with an even deeper stopband? Am I getting cumulative error and need double precision? Am I missing something obvious?

Since last night, I've made the following additional observations:

  • Changing everything double precision. including the original sinusoid generation np.exp(2j * np.pi * f * t) brought the stopband down below -90dB.
  • I woke up this morning thinking this might be the result of just the RMS calculation, as the numbers get very small after squaring, but converting to double just before the calculation doesn't help.
  • It also doesn't help to convert the original signal to double before applying the filter (which is also double).
  • With 23 bits of single precision, I must be able to do better than 60dB.

Here are the filter coefficients:

        coefficients = np.array([
                 0.0000060016235398,  0.0000225232064574, -0.0000090761946080, -0.0000185253128763,  0.0000054919335483,  0.0000280194533770,
                -0.0000026626363351, -0.0000385283331668, -0.0000034431756503,  0.0000504987635930,  0.0000135938348389, -0.0000632377133690,
                -0.0000287859523798,  0.0000757618568015,  0.0000499315016862, -0.0000867110007680, -0.0000777654415661,  0.0000943459911465,
                 0.0001127423719451, -0.0000965669022579, -0.0001549159433094,  0.0000909466345029,  0.0002038152171414, -0.0000748081392044,
                -0.0002583085688147,  0.0000453555855252,  0.0003164921863662,  0.0000001588307440, -0.0003756336075665, -0.0000642669244212,
                 0.0004321166167121,  0.0001490774184329, -0.0004814054413283, -0.0002559974896079,  0.0005181264333150,  0.0003854812580589,
                -0.0005362065594552, -0.0005368424812438,  0.0005289318003337,  0.0007078805812470, -0.0004893347868947, -0.0008948172631002,
                 0.0004103796416549,  0.0010920302321687, -0.0002854030332083, -0.0012919983382830,  0.0001084935034322,  0.0014852508773895,
                 0.0001250484221723, -0.0016604395232509, -0.0004182464428643,  0.0018045048046239,  0.0007719692587412, -0.0019029351357939,
                -0.0011844677730397,  0.0019401149807281,  0.0016509691466624, -0.0018997954472105, -0.0021633467846184,  0.0017656594972565,
                 0.0027098387951992, -0.0015219304042146, -0.0032748575978586,  0.0011540294011699,  0.0038389655782451, -0.0006492835473763,
                -0.0043789386395778, -0.0000023645806902,  0.0048678600770447,  0.0008077519284957, -0.0052754317309818, -0.0017698598225236,
                 0.0055682884655955,  0.0028873323308986, -0.0057102262479144, -0.0041541478810996,  0.0056625940124443,  0.0055593503835934,
                -0.0053842512098914, -0.0070870699304261,  0.0048314370898768,  0.0087166143844400, -0.0039568943452959, -0.0104228474662895,
                 0.0027081634880498,  0.0121766773521897, -0.0010242178022577, -0.0139457689059403, -0.0011706267819046,  0.0156953971820362,
                 0.0039777563460445, -0.0173894074278104, -0.0075447717571554,  0.0189912960389074,  0.0121054907843508, -0.0204653426305610,
                -0.0180639928239423,  0.0217777308244607,  0.0261908916184348, -0.0228976541697665, -0.0381558777897443,  0.0237983447290527,
                 0.0583374164208512, -0.0244579667873972, -0.1028686994496340,  0.0248603599090254,  0.3172249823495720,  0.4750043916617910,
                 0.3172249823495720,  0.0248603599090254, -0.1028686994496340, -0.0244579667873972,  0.0583374164208512,  0.0237983447290527,
                -0.0381558777897443, -0.0228976541697665,  0.0261908916184348,  0.0217777308244607, -0.0180639928239423, -0.0204653426305610,
                 0.0121054907843508,  0.0189912960389074, -0.0075447717571554, -0.0173894074278104,  0.0039777563460445,  0.0156953971820362,
                -0.0011706267819046, -0.0139457689059403, -0.0010242178022577,  0.0121766773521897,  0.0027081634880498, -0.0104228474662895,
                -0.0039568943452959,  0.0087166143844400,  0.0048314370898768, -0.0070870699304261, -0.0053842512098914,  0.0055593503835934,
                 0.0056625940124443, -0.0041541478810996, -0.0057102262479144,  0.0028873323308986,  0.0055682884655955, -0.0017698598225236,
                -0.0052754317309818,  0.0008077519284957,  0.0048678600770447, -0.0000023645806902, -0.0043789386395778, -0.0006492835473763,
                 0.0038389655782451,  0.0011540294011699, -0.0032748575978586, -0.0015219304042146,  0.0027098387951992,  0.0017656594972565,
                -0.0021633467846184, -0.0018997954472105,  0.0016509691466624,  0.0019401149807281, -0.0011844677730397, -0.0019029351357939,
                 0.0007719692587412,  0.0018045048046239, -0.0004182464428643, -0.0016604395232509,  0.0001250484221723,  0.0014852508773895,
                 0.0001084935034322, -0.0012919983382830, -0.0002854030332083,  0.0010920302321687,  0.0004103796416549, -0.0008948172631002,
                -0.0004893347868947,  0.0007078805812470,  0.0005289318003337, -0.0005368424812438, -0.0005362065594552,  0.0003854812580589,
                 0.0005181264333150, -0.0002559974896079, -0.0004814054413283,  0.0001490774184329,  0.0004321166167121, -0.0000642669244212,
                -0.0003756336075665,  0.0000001588307440,  0.0003164921863662,  0.0000453555855252, -0.0002583085688147, -0.0000748081392044,
                 0.0002038152171414,  0.0000909466345029, -0.0001549159433094, -0.0000965669022579,  0.0001127423719451,  0.0000943459911465,
                -0.0000777654415661, -0.0000867110007680,  0.0000499315016862,  0.0000757618568015, -0.0000287859523798, -0.0000632377133690,
                 0.0000135938348389,  0.0000504987635930, -0.0000034431756503, -0.0000385283331668, -0.0000026626363351,  0.0000280194533770,
                 0.0000054919335483, -0.0000185253128763, -0.0000090761946080,  0.0000225232064574,  0.0000060016235398,
            ]

and here is the code that applies them to the signal. Although my project is in C, this is in Python, as are all the plots, which I'm using to try things out. The output of the filter is the same, though.

return convolve(signal, coefficients, mode='same', method='fft')

I have used both numpy and scipy convolution. The above one, with the method='fft' is the scipy version. I was hoping the fft convolution would reduce cumulative error. It may have done so but it has no visible effect on the result.

I have set up the Python program to be able to select either single or double precision floating point. My intention was to use single precision FP in the C version. I'm using ordinary CPUs like X86_64 and ARM so floating point is comparable to integer performance. Double precision, though, would be a performance penalty I'd rather avoid.

I do not use rounding, truncation or any specific accumulator.

$\endgroup$
12
  • $\begingroup$ Could you clarify what are the units on the X axis in the frequency plots? $\endgroup$
    – MBaz
    Feb 1 at 23:10
  • 1
    $\begingroup$ x-axis for the filter is relative to the coefficients of the filter (215) and the FFT was zero padded to 512 points, the other FFTs are in kHz, and for the RMS plot in seconds. All plots y-axis are dB, except the last one. The original signal was +/- 0.9, so just under 0 dB. $\endgroup$
    – Deepstop
    Feb 1 at 23:15
  • 1
    $\begingroup$ Done. I've also started looking at double precision. A closer look at the original single precision FFT shows the noise about 100 dB down. I wonder if cumulative error is creeping in, although that's a lot of cumulative error. $\endgroup$
    – Deepstop
    Feb 2 at 3:43
  • 2
    $\begingroup$ It would help if you post details of your actual filter coefficients and the code to implement the filter. Do you use a double precisions accumulator? How often in the algorithm do you round or truncate? What exact rounding algorithm do you use ? Fixed point processing is all about the details $\endgroup$
    – Hilmar
    Feb 2 at 14:11
  • 1
    $\begingroup$ I ran your filter and what I assume your test signal is and got perfectly fine results. The wave form is perfectly smooth with small spikes at the frequency boundaries (which could be mitigated by ensuring phase continuity). I think something is wrong in your process, $\endgroup$
    – Hilmar
    Feb 2 at 15:35

2 Answers 2

2
$\begingroup$

I was trying to reproduce, but didn't see the errors.

Let's take a look at the filtered waveform

enter image description here

We see the shape as expected but there are bunch of spikes in there. These spikes are due to the discontinuity of the frequency. The hard switch triggers the transient responses of the filter.

If we calculate the the RMS in chunks of 2000 we get this

enter image description here

So get about 90 dB of attenuation EXCEPT in the chunks that contain a transient spike. This pattern depends on how the RMS chunks line up with the frequency changes.

An easy way around this would be to use a linear sweep instead of discrete frequencies. In this case, the waveform will have no spikes and the RMS picture looks like expected.

enter image description here

In this case both input and output were quantized to 14-bits but the filter itself and the convolution were done in double precision floating point (simply because that's the default data type in MATLAB).

$\endgroup$
2
$\begingroup$

I'm still not 100% clear but have a good idea of what's going on. I discovered that the problem lies with the test signals. These are on the following frequencies in the -16kHz to +16kHz range:

[-14750.0, -13750.0, -12750.0, -11750.0, -10750.0, -9750.0, -8750.0, 
  -7750.0,  -6750.0,  -5750.0,  -4750.0,  -3750.0, -2750.0, -1750.0, 
   -750.0,    250.0,   1250.0,   2250.0,   3250.0,  4250.0,  5250.0, 
   6250.0,   7250.0,   8250.0,   9250.0,  10250.0, 11250.0, 12250.0, 
  13250.0,  14250.0,  15250.0]

In the tests in the OP, each frequency was generated for 1/2 second at a time and written to a file. I noticed today that if I generated the test signal using double precision, even though the filter, FFT & RMS calculations were still single precision, I would get the >90dB stopband attenuation I was expecting. In fact, I would still get >90dB if I converted the signal to single precision after it was generated.

To get >90dB, I needed to use double precision for both the time intervals and the signal itself (exp(2j * pi * f * t).

I reduced the code to create the test signal to the following:

channel_spacing = 1000
offset = 250
carrier_amplitude = .9
pos = np.arange(0, (sample_rate - channel_spacing) / 2, channel_spacing, dtype=np.float64)
channels = np.concatenate((np.flip(-pos[1:]), pos )) +  offset

num_samples = int(dwell_seconds * channels.size * sample_rate)
t =  np.linspace(0, dwell_seconds * channels.size - 1, num_samples, dtype=np.float64)

wavec  = carrier_amplitude * np.exp(2j * np.pi * 
    channels.reshape(channels.size, -1) * t.reshape(channels.size, -1), 
    dtype=np.complex128).flatten().astype(np.complex64)

The weird thing is that the errors are cumulative. In the OP, each segment showed step decreases in stopband attenuation through any 1/2 second segment. In this simplified code, which uses a single array for t, these decreases in attenuation continue through the entire set of segments, reducing to about -30dB at the end. I've shown that the only difference between each value of t is only due to the difference in the number of significant digits, and I don't see a place where errors can accumulate.

To illustrate the last point, here's a test I did with t, which outputs zero.

t =   np.linspace(0, dwell_seconds * channels.size - 1, num_samples, dtype=np.float64)
tf =  np.linspace(0, dwell_seconds * channels.size - 1, num_samples, dtype=np.float32)
print(np.max(np.abs(t.astype(np.float32) - tf)))

I will probably take this over to Stack Overflow as it seems to be more of a programming question than DSP, but I thought I'd report back on my progress as a (partial) answer, rather than adding this to an already long question. I'll leave it unchecked in case anyone can take the answer a bit further than this.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.