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? _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
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 [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
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.----- 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 <[hidden email]> wrote:
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
Are these the proper frequency domain plots+ w2, h2 = freqz(data) 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()
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <[hidden email]> wrote:
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
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user freqresp.png (64K) Download Attachment |
Any idea why this plot shows positive gain?
Like +50dB around 500Hz? -Matti
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
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) 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 which doesn’t look like a highpass?
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
On Sun, Jul 10, 2016 at 12:00 PM, Matti Viljamaa <[hidden email]> wrote:
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
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
But what if I want to plot the effect of the lpf on the signal, rather than the filter magnitude response? -Matti
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
On Sun, Jul 10, 2016 at 12:08 PM, Matti Viljamaa <[hidden email]> wrote:
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) ----------
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user figure_3.png (38K) Download Attachment |
On Sun, Jul 10, 2016 at 12:57 PM, Warren Weckesser <[hidden email]> wrote:
Also, matplotlib has its own tools for this. See, for example, the `psd` function demonstrated in this example: http://matplotlib.org/examples/pylab_examples/psd_demo.html Warren
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Warren Weckesser-2
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] lpf = remez(513, bands, desired) gives the following plots: 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 [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
On Mon, Jul 11, 2016 at 7:03 AM, Matti Viljamaa <[hidden email]> wrote:
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
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
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.] -Matti
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
What about the following:
is this (the almost -15dB cut in the transition from <2000Hz to >2000Hz) a “feature" of the filter? plotted with: 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 desired = [0.5, 1, 0.3] -Matti _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Matti Viljamaa
On Mon, Jul 11, 2016 at 10:13 AM, Matti Viljamaa <[hidden email]> wrote:
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
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Matti Viljamaa
On Mon, Jul 11, 2016 at 10:23 AM, Matti Viljamaa <[hidden email]> wrote:
That dip in the gain looks like it is in the "don't care" interval [cutoff - 0.5*transition_width, cutoff + 0.5*transition_width] (i.e. [1900, 2100]). The remez algorithm doesn't specify the behavior in the "don't care" intervals. If it turns out that you get something desirable, then sure, it's a feature, but if you tweak the parameters (i.e. the band edges or filter length), the behavior in the "don't care" intervals could be completely different. Warren
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Matti Viljamaa
Just a suggestion: for low-pass and band-pass filters, you might have better luck with Butterworth filters:
from scipy.signal import butter, freqz fs = 44100 nyquist = 0.5/fs passband = np.array([2000.,4500.]) * nyquist order = 4 num, den = butter(order, passband, btype='bandpass') w, h = freqz(num, den) -clancy > On Jul 11, 2016, at 10:23 AM, Matti Viljamaa <[hidden email]> wrote: > > What about the following: > > is this (the almost -15dB cut in the transition from <2000Hz to >2000Hz) a “feature" of the filter? > > <Screen Shot 2016-07-11 at 17.22.02.png> > > plotted with: > > 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 > desired = [0.5, 1, 0.3] > > -Matti > _______________________________________________ > SciPy-User mailing list > [hidden email] > https://mail.scipy.org/mailman/listinfo/scipy-user _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Warren Weckesser-2
Are there general ways for ensuring that the remez produces a filter that’s “correct”, i.e. no huge undesired boosts or other such “glitches”?
-Matti _______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
On Tue, Jul 12, 2016 at 10:42 AM, Matti Viljamaa <[hidden email]> wrote: Are there general ways for ensuring that the remez produces a filter that’s “correct”, i.e. no huge undesired boosts or other such “glitches”? Not that I know of. I always check the filter response with freqz. Warren
_______________________________________________ SciPy-User mailing list [hidden email] https://mail.scipy.org/mailman/listinfo/scipy-user |
Free forum by Nabble | Edit this page |