# Trying to apply Parks-McClellan optimal FIR filter design

My goal is to apply the Remez exchange algorithm in R to design an FIR band pass filter for a specific frequency window. The question is not supposed to be referring to R of course but about how to choose the parameters correctly and why the result looks like as it does.

The toy signal is an additive combinations of four sine waves each adding a single frequency component to the final signal. The frequencies are 1/20, 1/300, 1/5000 and 1/7000. Its FFT plot is shown in figure (s).

Now what I apparently didn't get yet is how I would have to choose the parameters for the filter design - being the pass band and the stop band cutoff points.

The 1st filter uses the cutoff points:

• 0
• 1/10 - 2/1000
• 1/10 - 1/1000
• 1/10 + 1/1000
• 1/10 + 2/1000
• 1

With 1 being the Nyquist-frequency and the valid cutoff point interval being (0,1).

As you can see in figure (r1) this results in an attenuation of all frequency components except for 1/5000. (I would assume that also the frequency component 1/5000 is actually attenuated but relatively less then the other components)

So my first set of questions would be - why is that so? How could I calculate f.x. the frequency located in the center of the pass band? What is the quantitative relation between cutoff points on an interval (0,1) and the pass band?

The 2nd filter uses the cutoff points:

• 0
• 1/10 - 2/500
• 1/10 - 1/500
• 1/10 + 1/500
• 1/10 + 2/500
• 1

And the result is an (relative) attenuation of all components except for 1/300. What I am wondering about in this case is why the pass band is apparently completely different from the above filter even though I just made the filter's pass band wider. I would assume at least that 1/5000 is still located within the pass band and not attenuated. A possible explanation might be that 1/300 is now included within the pass band and relative to 1/5000 less attenuated making 1/5000 look just smaller.

I would like to add that I am - as you can see - quite new to this subject and actually from a different area and trying my best to understand what's going on.

The canonical question would be: How do I have to choose the cutoff points in case of a possible interval (0,1=Nyquist-Freq) to attenuate all frequency components except for 1/N and its close neighbourhood?

library(signal)
par(mfrow=c(5,1))

# the toy signal
ch <- sin(2*pi*1:100000/floor(100000/20)) +
sin(2*pi*1:100000/floor(100000/300)) +
sin(2*pi*1:100000/floor(100000/5000)) +
sin(2*pi*1:100000/floor(100000/7000))

# (s)
barplot(abs(fft(ch))[1:10000])

filter <- function(c0,d1,d2) {

# the filter design with cutoffs specified
fir <- remez(n=1000,
f=c(0,c0-d2,c0-d1,c0+d1,c0+d2,1),
a=c(0,0,1,1,0,0))

# (fN)
plot(freq$f,abs(freq$h),type="l")
freq <- freqz(fir,n=100000)
y <- signal::filter(as.vector(fir), 1, x=ch)

# (sN)
barplot(abs(fft(y))[1:10000])
}

filter(1/10,1/1000,2/1000)
filter(1/10,1/500,2/500)


## 1 Answer

There are a few problems with your designs. First, the desired bandwidths and transition bands (between passbands and stopbands) are very narrow. Second, some of the center frequencies are very low, especially the ones at 20/100000 and 300/100000. And finally, your chosen filter length is very large. I do not know the implementation of the Remez exchange algorithm that you use, but I've seen several implementations that run into serious numerical problems for such filter lengths. You can see in your designs that around Nyquist there are some small peaks in the filter response. These are due to numerical problems.

So you should try to make the passband and the transition bands a bit wider and reduce the filter length to see if you can get a reasonable design. From there you could try to make the specifications a bit more stringent and increase the filter order and see how far you can go.

In Octave, the following specification gives a reasonable result (center frequency: 0.1, filter length: 300, '1' corresponds to Nyquist):

h=remez(300,[0,.07,.085,.115,.13,1],[0,0,1,1,0,0]);

• Thank you very much for your answer! Does that mean that my approach is at least theoretically correct? Could you give me a (theoretically) correct formula on how to calculate the cutoffs? Currently my idea is that the center of the cutoffs is essentially the frequency i.e. 1/7000 for a pass band aiming at 7000 periods per full signal (one might say 7000 Hz if the signal takes a second) and the cutoffs are supposed to be chosen very close to that middle point - ignoring potential numerical issues. – Raffael May 13 '14 at 17:32
• There is one point I especially don't get. If I f.x. want to let frequencies around 1000 Hz pass, then the center between the two pass band cutoffs would be - if my idea is correct - 1/1000. Then the transition bands are inevitably very small f.x. Because I would have to fit a stop and a pass band cutoff between 0 and 1/1000. That is one thought that leads me to thinking there is something wrong with my understanding of how this works. – Raffael May 13 '14 at 17:44
• @Raffael: Well, you know that in the discrete-time domain you work with relative frequencies. All frequencies are normalized by the sampling frequency. Your frequency '7000' is actually 2*7000/100000 (if Nyquist is '1'). A sinusoidal signal is $x(n)=\sin(2\pi nf/f_s)$ where $f$ is the frequency in Hertz, and $f_s$ is the sampling frequency in Hertz. If you compare this to your signal, you can see the relation between your numbers and the relative frequency. – Matt L. May 13 '14 at 17:48
• NICE! Now I get what I was looking for! Thank you very much – Raffael May 13 '14 at 17:58