🚀 go-pugleaf

Howdy,

Yep, I shamefully accept the brain dead award of comp.dsp for as Miguel  
points out:
> 
> I believe the problem is not in the filter design, but in the
> generation of the averaged noise spectra.  I noticed in the Matlab
> script included in your posting that you did not window the filtered 
> data in the calculation of the averaged spectra.  
> This results in unacceptably high "spectral leakage" into the
> stopband region....

the filter is operating properly, it is just the cyclic nature of the the  
FFT (used to examing the resulting spectrum) which dislikes the  
discontinuity produced when y(N-1) wraps around to y(0). From the FFT  
perspective, such a discontinuity requires lots of high frequency  
components. By applying a window to the data set, we reduce the energy  
associated with this discontinuity and the FFT-based spectrum converges on  
the true (infinite time) spectrum. I believe a similar convergence would  
also occur if I used very long data sets so that the spectral energy  
associated with this discontinuity was reduced below the noise floor of  
the filter. Any comments?

Additionally, John Fleming and Tony Richardson found an error in the  
Matlab script; the line that reads,

    [y,zf]=filter(remezCoeff,1,x,yf);

should be

    [y,yf]=filter(remezCoeff,1,x,yf);

so that if we combine this correction with Miguel Oyarzun's recommendation  
we have:

> The following changes to your script will illustrate the problem.
> 
> >% 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);
> 
>    Wind = blackman(N);              
> 
> >  for index=1:DataBlocks
> >    x=rand(1,N);   % generate the data set of N points
> >    [y,yf]=filter(remezCoeff,1,x,yf);   % filter the data (error here)
> >    [z,zf]=filter(hammingCoeff,1,x,zf);
> 
>      Rspectrum = Rspectrum + abs(fft(Wind'.*y)); % add the results
>                                %     ^^^^
>      Hspectrum = Hspectrum + abs(fft(Wind'.*z));
>                                %     ^^^^
> >  end
> 

Using this modification, the averaged spectra from pumping white noise  
into the designed filter does, in fact, closely mimic that of the designed  
filter!!

Thanks to Miguel, John, and Tony for their efforts.

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