implementaion of remez exchange algorithm for filter design

Viewed 822

I have found C implementation for remez exchange algorithm(https://github.com/janovetz/remez-exchange/blob/master/remez.c) and I want to use this to implement a generic prototype lowpass filter. Although the comments have been provided to guide on how to use it I am having trouble implementing it. Could anyone tell me what are the correct inputs which I am suppose to give to the remez() function to get coefficients for lowpass filter because I suspect I am not providing the correct inputs.

Some information regarding what I am trying to do: I have been reading and experimenting with the code and I found out that for lowpass filter I should give type as 1(i.e Bandpass). Anyhow just to try the code I gave the parameters for bandpass filter in MATLAB and same parameters to remez() but the results are very different. Could somebody have a look at my coefficients and tell me what I am doing wrong because I don't get it?

MATLAB:

% All frequency values are in Hz.
Fs = 2000;  % Sampling Frequency
Fstop1 = 200;              % First Stopband Frequency
Fpass1 = 300;              % First Passband Frequency
Fpass2 = 700;              % Second Passband Frequency
Fstop2 = 800;              % Second Stopband Frequency
Dstop1 = 0.000177827941;   % First Stopband Attenuation
Dpass  = 0.0057563991496;  % Passband Ripple
Dstop2 = 0.000177827941;   % Second Stopband Attenuation
dens   = 20;               % Density Factor
% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fstop1 Fpass1 Fpass2 Fstop2]/(Fs/2), [0 1 ...
                           0], [Dstop1 Dpass Dstop2]);

These values (N, Fo, Ao, W) obtained from MATLAB, I feed to remez() function.

So the C++ function looks like: (just the part where I am calling the remez function)

int numtaps = 64;
int numbands = 3;
std::vector<double> h(numtaps);
double  bands[] = [0, 0.2000, 0.3000, 0.7000, 0.8000, 1.000];
double  des[] = [0, 0, 1, 1, 0, 0];
double weight[] = [0.000177827941, 0.0057563991496, 0.000177827941];
int type = 1;
int griddensity = 20;
remez(&(h[0]), int *numtaps, int *numband, const double bands[], const double des[], const double weight[], int *type, int *griddensity);

Could anyone please tell me how to use this code properly. I have not been able to find anything online about it.

1 Answers

Matlab's firpm and firpmord use band edge frequencies normalized by the Nyquist rate, such that the specified frequencies fall in the [0,1] range. The remez function implementation you used on the other hand assumes that band edge frequencies are normalized by the sampling rate, such these frequencies fall in the [0,0.5] range. As a result you should divide each specified frequency by 2:

double  bands[] = [0, 0.2000/2, 0.3000/2, 0.7000/2, 0.8000/2, 1.000/2];

Also, to have smaller ripples in a given band you must increase the weight associated with that band, so to get the correct relative ripple amplitudes you should provide weights which are inversely proportional to the ripples:

double weight[] = {1.0/0.000177827941, 1.0/0.0057563991496, 1.0/0.000177827941};

You should then be able to call remez as follow:

remez(&(h[0]), &numtaps, &numbands, bands, des, weight, &type, &griddensity);

This should provide you a reasonable approximation of the specified bandpass filter. Note that since the desired filter characteristics are provided in the frequency-domain, you should be comparing the frequency responses of the designed filter rather than the filter coefficients h (which are in the time-domain). You could freqz for this purpose. To actually meet the filter specification, you will probably have to increase N from the estimated N=64 given by firpmord to around 68.

Finally, since you mentioned that your goal is to obtain a lowpass filter, you would achieve this by specifying a bandpass filter with 2 bands, with the first band having a desired value of 1 and the second a desired value of 0, in a way similar to the following:

int numbands = 2;
double  bands[] = [0, 0.2000/2, 0.3000/2, 0.5];
double  des[]   = [1, 1, 0, 0];
...
Related