[SciPy-User] scipy.signal.remez producing odd filter effects / not producing desired effect

Mon Jul 11 11:50:57 EDT 2016

On Mon, Jul 11, 2016 at 10:13 AM, Matti Viljamaa <mviljamaa at kapsi.fi> wrote:

> On 11 Jul 2016, at 16:30, Warren Weckesser <warren.weckesser at gmail.com>
> wrote:
>
>
>
> On Mon, Jul 11, 2016 at 7:03 AM, Matti Viljamaa <mviljamaa at kapsi.fi>
> wrote:
>
>>
>> On 10 Jul 2016, at 19:57, Warren Weckesser <warren.weckesser at gmail.com>
>> wrote:
>>
>>
>>
>> On Sun, Jul 10, 2016 at 12:08 PM, Matti Viljamaa <mviljamaa at kapsi.fi>
>> wrote:
>>
>>>
>>> But what if I want to plot the effect of the lpf on the signal, rather
>>> than the filter magnitude response?
>>>
>>>
>>
>> For that, you want to plot the spectral content of the original and
>> filtered data.  You can use a Fourier transform and plot the magnitude of
>> the Fourier coefficients against the frequencies, or you can use a function
>> such as `scipy.signal.periodogram` or `scipy.signal.welch` which take care
>> of the FFT details for you.
>>
>> Here's a new version of my script.  It now creates a third figure
>> containing plots of the periodogram of the original data and the filtered
>> data.  I used `scipy.signal.periodogram`, but I recommend experimenting
>> with `welch` also.
>>
>> Warren
>>
>> ----------
>>
>> from __future__ import division
>>
>> from scipy.signal import remez, freqz, lfilter, periodogram
>> from scipy.io import wavfile
>> import numpy as np
>> import matplotlib.pyplot as plt
>>
>>
>> fs = 44100
>>
>> # Design a low-pass filter using remez.
>> cutoff = 2000.0
>> transition_width = 200
>> bands = np.array([0, cutoff - 0.5*transition_width,
>>                   cutoff + 0.5*transition_width, fs/2.0]) / fs
>> desired = [1, 0]
>> lpf = remez(513, bands, desired)
>>
>> # Plot the frequency response of the filter.
>> w, h = freqz(lpf)
>> plt.figure(1)
>> plt.plot(fs*w/(2*np.pi), 20*np.log10(abs(h)))
>> plt.xlim(0, fs/2)
>> plt.xlabel('Frequency (Hz)')
>> plt.ylabel('Gain (dB)')
>> plt.grid(True)
>>
>> # Create a test signal with three frequencies: two inside the pass-band,
>> # are one far outside the pass-band that should be filtered out.
>> T = 0.5
>> nsamples = int(T*fs)
>> t = np.linspace(0, T, nsamples, endpoint=False)
>> freq = 440
>> data = np.sin(2*np.pi*freq*t) + 0.5*np.cos(2*np.pi*(2*freq)*t)
>> data += 0.75*np.sin(2*np.pi*(cutoff + 5*transition_width)*t)
>> data /= 1.01*np.abs(data).max()
>>
>> # Filter the input using lfilter. (Alternatively, convolution could be
>> used.)
>> filtered_data = lfilter(lpf, 1, data)
>>
>> # Plot the input and output in the same figure.
>> plt.figure(2)
>> plt.plot(t, data, 'b', label='original')
>> plt.plot(t, filtered_data, 'g', label='filtered')
>> plt.xlabel('Time (sec)')
>> plt.ylim(-1.1, 1.1)
>> plt.legend()
>>
>> plt.figure(3)
>> freqs, data_spec = periodogram(data, fs=fs)
>> freqs, filtered_data_spec = periodogram(filtered_data, fs=fs)
>> plt.subplot(2, 1, 1)
>> plt.plot(freqs, data_spec, 'b', label='original')
>> plt.xlim(0, 3*cutoff)
>> plt.axvline(cutoff, color='k', alpha=0.15)
>> plt.legend()
>> plt.subplot(2, 1, 2)
>> plt.plot(freqs, filtered_data_spec, 'g', label='filtered')
>> plt.xlim(0, 3*cutoff)
>> plt.axvline(cutoff, color='k', alpha=0.15)
>> plt.legend()
>> plt.xlabel('Frequency (Hz)')
>>
>> plt.show()
>>
>> # Save the test signal and the filtered signal as wav files.
>> wavfile.write("data.wav", fs, data)
>> wavfile.write("filtered_data.wav", fs, filtered_data)
>>
>>
>> —————
>>
>>
>> Thanks for this. However, I think there’s some other problem.
>>
>> Modifying the above remez to the following:
>>
>> fs = 44100
>>
>> # Design a low-pass filter using remez.
>> cutoff = 2000.0
>> transition_width = 200
>> bands = np.array([0, cutoff - 0.5*transition_width,
>>                   cutoff + 0.5*transition_width, cutoff +
>> 0.5*transition_width + 1200.0, cutoff + 0.5*transition_width +
>> 2400.0,fs/2.0]) / fs
>> desired = [0, 1, 0]
>>
>
>
> You specified three intervals in `bands`: [0, 1900], [2100, 3300] and
> [4500, 22050] with gains 0, 1, and 0, respectively (i.e. you are designing
> a bandpass filter).  The remez algorithm designs a filter that is
> equiripple in the specified bands.  Outside those bands, however, the
> filter gain is unspecified; typically these intervals are described as the
> "don't care" intervals.  The remez algorithm makes no promises about the
> gain in these intervals, and it can--as you have found--result in
> unacceptable behavior.  The big peak in the filter gain occurs in the
> "don't care" interval [3300, 4500].
>
> To fix it, you can try reducing the length of that "don't care" interval,
> or increasing the length of the filter (if that is an option).  It might
> require several iterations of band adjustments to achieve a desirable
> result.  You can also try the `weight` argument--perhaps give the stop
> bands a lower weight than the pass band.  Some experimentation and
> iteration will be involved.
>
> Warren
>
>
> So doing e.g.
>
> bands = np.array([0, cutoff - 0.5*transition_width,
>                   cutoff + 0.5*transition_width, cutoff +
> 0.5*transition_width + 2400.0, cutoff + 0.5*transition_width +
> 2500.0,fs/2.0]) / fs
>
> i.e.
>
> [     0.,   1900.]   [2100.,   4500.]   [4600.,  22050.]
>
> the filter is shaped properly.
>
> However,
>
> do you know why
>
> [     0.,   1900.]   [2100.,   4500.]   [4500.,  22050.]
>
> fails? I.e. when the transition band is 0Hz?
>
>

You are asking for an ideal transition at 4500 Hz from perfect pass band to
perfect stop band.  The remez algorithm doesn't converge.  You'd have to
dig into the details of the algorithm for a good answer to "why", so I'll
just say that the problem is too hard--possibly even impossible--for the
remez algorithm.

Warren

> -Matti
>
>
> lpf = remez(513, bands, desired)
>>
>> gives the following plots:
>>
>> <Screen Shot 2016-07-11 at 14.01.40.png>
>>
>> Why do these parameters lead to that +150dB peak, the peak in the
>> “filtered" frequency plot and the strange full bandwidth burst in the last
>> picture?
>>
>> -Matti
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> https://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> https://mail.scipy.org/mailman/listinfo/scipy-user
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> https://mail.scipy.org/mailman/listinfo/scipy-user
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20160711/472249ee/attachment.html>