[SciPy-User] scipy.signal.remez producing odd filter effects / not producing desired effect
Warren Weckesser
warren.weckesser at gmail.com
Sun Jul 10 12:57:56 EDT 2016
On Sun, Jul 10, 2016 at 12:08 PM, Matti Viljamaa <mviljamaa at kapsi.fi> wrote:
>
> On 10 Jul 2016, at 19:04, Warren Weckesser <warren.weckesser at gmail.com>
> wrote:
>
>
>
> On Sun, Jul 10, 2016 at 12:00 PM, Matti Viljamaa <mviljamaa at kapsi.fi>
> wrote:
>
>> Or why doing
>>
>> 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 = [0, 1] # <------
>> lpf = remez(513, bands, desired)
>>
>> 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)
>> data += np.sin(2*np.pi*(cutoff + 5*transition_width)*t)
>> data /= 1.01*np.abs(data).max()
>>
>> filtered_data = lfilter(lpf, 1, data)
>>
>> w3, h3 = freqz(filtered_data)
>>
>>
>
> Matti,
>
> Did you see the email that I sent on July 6, and the plot that I attached
> to it? Here's what I said:
>
> The arguments to 'freqz' are the filter coefficients, not the signal. See
> the example script in my previous email, where I have
>
> w, h = freqz(lpf)
>
> Warren
>
>
>
> 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)
----------
> -Matti
>
>
>
>> plt.figure(3)
>> plt.plot(fs*w3/(2*np.pi), 20*np.log10(abs(h3)))
>> plt.xlim(0, fs/2)
>> plt.xlabel('Frequency (Hz)')
>> plt.ylabel('Gain (dB)')
>> plt.grid(True)
>>
>> plt. show()
>>
>> gives
>>
>> <Screen Shot 2016-07-10 at 18.58.09.png>
>>
>> which doesn’t look like a highpass?
>>
>>
>> On 07 Jul 2016, at 13:08, Matti Viljamaa <mviljamaa at kapsi.fi> wrote:
>>
>> Any idea why this plot shows positive gain?
>> Like +50dB around 500Hz?
>>
>> -Matti
>>
>> On 06 Jul 2016, at 17:57, Warren Weckesser <warren.weckesser at gmail.com>
>> wrote:
>>
>>
>>
>> On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <mviljamaa at kapsi.fi>
>> wrote:
>>
>>> Are these the proper frequency domain plots+
>>>
>>> w2, h2 = freqz(data)
>>>
>>>
>>
>> The arguments to 'freqz' are the filter coefficients, not the signal.
>> See the example script in my previous email, where I have
>>
>> w, h = freqz(lpf)
>>
>> The attached plot is the frequency response that I get when I run that
>> script.
>>
>> Warren
>>
>>
>>
>>
>>> plt.figure(2)
>>> plt.plot(fs*w2/(2*np.pi), 20*np.log10(abs(h2)))
>>> plt.xlim(0, fs/2)
>>> plt.xlabel('Frequency (Hz)')
>>> plt.ylabel('Gain (dB)')
>>> plt.grid(True)
>>>
>>> w3, h3 = freqz(filtered_data)
>>>
>>> plt.figure(3)
>>> plt.plot(fs*w3/(2*np.pi), 20*np.log10(abs(h3)))
>>> plt.xlim(0, fs/2)
>>> plt.xlabel('Frequency (Hz)')
>>> plt.ylabel('Gain (dB)')
>>> plt.grid(True)
>>>
>>> plt. show()
>>>
>>> <Screen Shot 2016-07-06 at 14.53.08.png>
>>>
>>> On 04 Jul 2016, at 18:45, Warren Weckesser <warren.weckesser at gmail.com>
>>> wrote:
>>>
>>> Matti,
>>>
>>> I don't have your input file, so I can't reproduce your result exactly.
>>> For a basic demonstration, it is probably easier to create a signal in the
>>> code instead of requiring an input wav file. Here's an example that does
>>> that. The input signal (called 'data') is the sum of two sine waves, one
>>> at 440 Hz and one at 3000 Hz. The cutoff frequency for the low-pass filter
>>> is 2000 Hz. (Actually the transition from the pass-band to the stop-band is
>>> from 1900 Hz to 2100 Hz.) You should be able to see in the plot and hear
>>> in the wav files that the high frequency component has been filtered out.
>>>
>>> Warren
>>>
>>> -----
>>>
>>> from __future__ import division
>>>
>>> from scipy.signal import remez, freqz, lfilter
>>> 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 two frequencies: one 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)
>>> data += 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)
>>> plt.plot(t, filtered_data)
>>>
>>> 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)
>>>
>>>
>>> -----
>>>
>>>
>>>
>>>
>>> On Mon, Jul 4, 2016 at 9:52 AM, Matti Viljamaa <mviljamaa at kapsi.fi>
>>> wrote:
>>>
>>>>
>>>> On 04 Jul 2016, at 14:00, Matti Viljamaa <mviljamaa at kapsi.fi> wrote:
>>>>
>>>> I’m trying to use the scipy.signal.remez filter design function. What
>>>> I’m doing is the following (with the help of
>>>> http://www.ee.iitm.ac.in/~nitin/teaching/ee5480/firdesign.html):
>>>>
>>>> import os
>>>> import scipy.io.wavfile
>>>> from scipy import signal
>>>> from pylab import *
>>>>
>>>> os.chdir("/Users/mviljamaa/Music")
>>>> sr, data = scipy.io.wavfile.read(open(“sound.wav", 'r'))
>>>>
>>>> lpf = signal.remez(21, [0, 0.05, 0.1], [1.0, .0])
>>>>
>>>> # Plot the magnitude response
>>>> w, h = signal.freqz(lpf)
>>>> plot(w/(2*pi), 20*log10(abs(h)))
>>>> show()
>>>>
>>>> # Filtered data
>>>> sout = signal.lfilter(lpf, 1, data)
>>>>
>>>> scipy.io.wavfile.write("equalized.wav", sr, sout)
>>>>
>>>> —
>>>>
>>>> Now judging by the magnitude response of the filter lpf, it should be a
>>>> lowpass of some sort (I’m not sure how to interpret the cutoff frequency
>>>> yet).
>>>>
>>>> The output I get is substantially gained (up to the point where it
>>>> sounds distorted) and I cannot hear the lowpass filter.
>>>>
>>>> What am I doing wrong?
>>>>
>>>>
>>>>
>>>> Using a slightly modified code (from http://pastebin.com/LPEtXdzx):
>>>>
>>>> import os
>>>> import scipy.io.wavfile
>>>> from scipy import signal
>>>> from pylab import *
>>>> import numpy as np
>>>>
>>>> os.chdir("/Users/mviljamaa/Music")
>>>> sr, data = scipy.io.wavfile.read(open(“sound.wav", 'r'))
>>>>
>>>> fs = 44100
>>>>
>>>> bands = array([0,3500,4000,5500,6000,fs/2.0]) / fs
>>>> desired = [1, 0, 0]
>>>> lpf = signal.remez(513, bands, desired)
>>>>
>>>> from scipy.signal import freqz
>>>> w, h = signal.freqz(lpf)
>>>> plot(w/(2*pi), 20*log10(abs(h)))
>>>> show()
>>>>
>>>> sout = signal.lfilter(lpf, 1, data)
>>>>
>>>> sout /= 1.05 * max(abs(sout))
>>>>
>>>> scipy.io.wavfile.write("equalized.wav", sr, sout)
>>>>
>>>>
>>>> —
>>>>
>>>> This one is able to retain the gain, but I still don’t hear any lowpass
>>>> filtering.
>>>>
>>>> _______________________________________________
>>>> 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
>>>
>>>
>> <freqresp.png>_______________________________________________
>> 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
>
>
>
> _______________________________________________
> 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/20160710/23487863/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: figure_3.png
Type: image/png
Size: 28659 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20160710/23487863/attachment.png>
More information about the SciPy-User
mailing list