🚀 go-pugleaf

Hi folks,

I'm hoping somebody out there in netland can help me with a little problem  
associated with the design of a lowpass filter to be used as part of a  
decimating FIR filter. I used Matlab's "remez" function to design a filter  
having 0.5 dB ripple in the passband and greater than 80 dB down in the  
stopband.

Just for giggles, I decided to generate a white noise data set and run it  
through the filter using a time-domain convolution. Although the resulting  
spectrum was as designed, this test wasn't really representative of the  
the actual case which would be processing a continuous data stream on a TI  
C30. Therefore, to avoid the "ramping" behavior introduced by the  
convolution, I just grabbed the interior section (i.e., after the filter  
had been loaded) and did an fft. ALAS, this spectrum did not match my  
filter design--in fact, it missed my criteria for stopband rejection by 3-  
or 40 dB!!!

If you use the following Matlab script, you'll see what the problem is.  
The script should display four plots. The first is the diplay of the  
designed filter response, the second is the spectrum due to averaging  
white noise sent through a convolution with the filter and taking the DFT.  
The third plot shows the results of continuously filtering the signal so  
that the convolution "ramp up and down" region isn't included in the DFT.  
(This uses the Matlab "filter" function rather than the original  
approach.) Note that although this is what will be happening in real-life,  
that the resulting spectra doesn't match what was designed!! In fact, the  
anti-aliasing performance of a Hamming window-based filter appears much  
better. The final plot simply shows the passband portion of plot 3 with  
the passband edge drawn in. (The "DataBlocks" parameters determines how  
many data sets are averaged--so that some detail into the filter  
performance may be garnered.)

I suspect that the problem with the filter stems from a filter design  
presumption of a cyclic data set; however, I can't prove it. Alternate  
hypotheses include that the continuous filtered spectra represent an  
convolution of a sinc function (due to a rectangular data set) with the  
designed frequency response. If you have any ideas as to what creates this  
behavior, I'd appreciate the edification. More importantly, I still need  
to design a filter having the desired behavior when operating on a  
continuous data stream!!!!

I've run this phenomena past a few folk here at Penn State and they  
haven't pointed out any obvious areas where I'm out to lunch--may one of  
you will do better. ;) Is this sort of behavior a feature of the  
remez-exchange algorithm, a feature of my implementation, a feature of  
Matlab, or just the product of a confused and dangerous mind? If anybody  
out there can help me, I'd appreciate it.

Thanks,

Mark.
---
Mark Kotanchek
Guidance & Control Dept - 363 ASB
Applied Research Lab/Penn State
P.O. Box 30
State College, PA 16804

e-mail:	mek@guinan.psu.edu (NeXTmail)
TEL:	(814)863-0682
FAX:	(814)863-7843

----------------------Matlab Script Starts Here-------------------------


% Matlab script to illustrate discrepency with the remez-exchange
% filter behavior with respect to behavior exhibited when white
% noise is run through the filter continuously. The problem is
% to design a low-pass filter appropriate for a decimating FIR 
% filter which avoids aliasing upon decimation. The control 
% parameters listed below are commented; however, "N" represents
% the size of a data block of used in doing the FFT while 
% "DataBlocks" is the number of data sets averaged to produce some
% of the plots.
%
% This script was written using Matlab 3.5 on a Macintosh Powerbook
% by Mark Kotanchek at ARL/Penn State (mek@guinan.psu.edu).
%
% In the plots, the remez filter behavior is illustrated by a solid
% line while for comparison, a Hamming window-based filter behavior
% is shown via a dashed line.

% Control Parameters specifying the processing
  N=128;            % number of samples in a data set
  n = 31;           % Number of Filter coefficients minus one!!
  D = 5;	    % Decimation factor
  PassBandWt = 1;   % Parameters for the remez exchange filter
  StopBandWt = 800; % design function.
  NyquistFreq = 28; % Folding Frequency (half the sampling rate)
  PassFreq = 3.5;   % Passband edge
  DataBlocks = 20;  % Number of data sets to average
  rand('normal');   % Noise will be N(0,1) random variables

% This section of the script calculates the normalized filter 
% frequencies such that 1.0 corresponds to half the sampling
% rate.
  StopFreq =PassFreq+2*(NyquistFreq/D - PassFreq); % Foldback Freq.
  fEdge = PassFreq/NyquistFreq; % fractional passband edge.
  fStop = StopFreq/NyquistFreq; % fractional stopband edge.

% Calculate the remez exchange filter coefficients
  f=[0 fEdge fStop 1];
  a=[1 1 0 0];
  wtx=[PassBandWt StopBandWt];
  remezCoeff=remez(n,f,a,wtx); % the designed filter's coefficients

% Calculate the window-based filter coefficients for comparison.
  hammingCoeff = fir1(n,fEdge);

% Calculate the frequency response of the specified filters
  N1=256;  % resolution of frequency response
  [remezResponse,freq] = freqz(remezCoeff,1,N1);
  remezResponse = 20.*log10(abs(remezResponse));
  hammingResponse = 20.*log10(abs(freqz(hammingCoeff,1,N1)));
  responseFreqs = freq*NyquistFreq/pi;

% Plot the two frequency responses
  clg;
  subplot(221);
 plot(responseFreqs,remezResponse,'-',responseFreqs,hammingResponse,'--');
  title(['Filter Responses for ',int2str(n+1),' filter coefficients']);
  xlabel('frequency (kHz)');
  ylabel('Power (dB)');
  grid;

% Compute the behavior based upon filter and data set convolution

  % Initialize the data vectors
    x=rand(1,N);   % generate the data set

    yRemez=conv(x,remezCoeff);   % convolve with Remez coefficients
    YRemez=abs(fft(yRemez));
    YRemez(length(YRemez)+1)=YRemez(1);
    convRemez=YRemez;

    yHamming=conv(x,hammingCoeff);  % ditto with Hamming coefficients
    YHamming=abs(fft(yHamming));
    YHamming(length(YHamming)+1)=YHamming(1);
    convHamming=YHamming;

  for index=2:DataBlocks
    x=rand(1,N);   % generate the data set

    yRemez=conv(x,remezCoeff);   % convolve with Remez coefficients
    YRemez=abs(fft(yRemez));
    YRemez(length(YRemez)+1)=YRemez(1);
    convRemez=convRemez+YRemez;

    yHamming=conv(x,hammingCoeff);  % ditto with Hamming coefficients
    YHamming=abs(fft(yHamming));
    YHamming(length(YHamming)+1)=YHamming(1);
    convHamming=convHamming+YHamming;
  end
  convFreqs=(0:length(yRemez))*(2*NyquistFreq)/length(yRemez);
  convRemez=20.*log10(convRemez./DataBlocks);
  convHamming=20.*log10(convHamming./DataBlocks);
  
  subplot(222);
  plot(convFreqs,convRemez,'-',convFreqs,convHamming,'--');
  title('convolved Filter spectra');

% Filter the signal against a real data stream
  x=rand(1,N);
  [y,yf]=filter(remezCoeff,1,x);   % load up the filter
  [z,zf]=filter(hammingCoeff,1,x); % load up the filter

  Rspectrum=zeros(1,N);
  Hspectrum=zeros(1,N);

  for index=1:DataBlocks
    x=rand(1,N);   % generate the data set of N points

    [y,zf]=filter(remezCoeff,1,x,yf);   % filter the data
    [z,zf]=filter(hammingCoeff,1,x,zf);

    Rspectrum = Rspectrum + abs(fft(y)); % add the results
    Hspectrum = Hspectrum + abs(fft(z));
  end

  % average the spectra and compute the power
    Rspectrum = 20.*log10(Rspectrum./DataBlocks);
    Rspectrum(N+1)=Rspectrum(1);
    Hspectrum = 20.*log10(Hspectrum./DataBlocks);
    Hspectrum(N+1)=Hspectrum(1);
    frequencies = (0:N)*(2*NyquistFreq)/N;

  % plot the averaged spectra
    subplot(223);
    plot(frequencies,Rspectrum,'-',frequencies,Hspectrum,'--');
    title('continuous Remez (solid)/Hamming (dashed)');
    xlabel('frequency (kHz)');
    ylabel('Power (dB)');
    grid;

  % plot some detail in the passband region
    stopIndex=round(fStop*(N+1));  % the stopband start index
    detailFreqs=frequencies(1:stopIndex);
    hammingDetail=Hspectrum(1:stopIndex);
    remezDetail=Rspectrum(1:stopIndex);
    subplot(224);
    plot(detailFreqs,remezDetail,'-',detailFreqs,hammingDetail,'--');
    limits=axis;
    hold on;
    passbandEdge=[PassFreq;PassFreq];
    passbandlines=[limits(3);limits(4)];
    plot(passbandEdge,passbandlines,'-.');
    hold off;
    title('behavior in passband region');
    xlabel('frequency (kHz)');
    ylabel('Power (dB)');
    grid;