python filter design__python

來源:互聯網
上載者:User











import numpy as np;import matplotlib.pyplot as pltimport cv2;import math;from cmath import sinfrom numpy.core.numeric import dtypefrom scipy.signal import butter, lfilter, freqzfrom scipy.fftpack import fft;'''    A few comments:The Nyquist frequency is half the sampling rate.You are working with regularly sampled data, so you want a digital filter, not an analog filter. This means you should not use analog=True in the call to butter, and you should use scipy.signal.freqz (not freqs) to generate the frequency response.One goal of those short utility functions is to allow you to leave all your frequencies expressed in Hz. You shouldn't have to convert to rad/sec. As long as you express your frequencies with consistent units, the scaling in the utility functions takes care of the normalization for you.Here's my modified version of your script, followed by the plot that it generates.'''def butter_lowpass(cutoff, fs, order=5):    nyq = 0.5 * fs    normal_cutoff = cutoff / nyq    b, a = butter(order, normal_cutoff, btype='low', analog=False)    return b, adef butter_lowpass_filter(data, cutoff, fs, order=3):    b, a = butter_lowpass(cutoff, fs, order=order)    y = lfilter(b, a, data)    return y# Filter requirements.order = 3fs = 30.0       # sample rate, Hzcutoff = 3.667  # desired cutoff frequency of the filter, Hz# Get the filter coefficients so we can check its frequency response.b, a = butter_lowpass(cutoff, fs, order)print b;print a;# Plot the frequency response.w, h = freqz(b, a, worN=8000)plt.figure(1);plt.subplot(2, 1, 1)plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b')#plt.plot(w, np.abs(h), 'b')plt.plot(cutoff, 0.5*np.sqrt(2), 'ko')plt.axvline(cutoff, color='k')plt.xlim(0, 0.5*fs)plt.title("Lowpass Filter Frequency Response")plt.xlabel('Frequency [Hz]')plt.grid()# Demonstrate the use of the filter.# First make some data to be filtered.T = 5.0         # secondsn = int(T * fs) # total number of samplest = np.linspace(0, T, n, endpoint=False)# "Noisy" data.  We want to recover the 1.2 Hz signal from this.data = np.sin(1.2*2*np.pi*t) + 1.5*np.cos(9*2*np.pi*t) + 0.5*np.sin(12.0*2*np.pi*t)# Filter the data, and plot both the original and filtered signals.y = butter_lowpass_filter(data, cutoff, fs, order)plt.subplot(2, 1, 2)plt.plot(t, data, 'b-', label='data')plt.plot(t, y, 'g-', linewidth=2, label='filtered data')plt.xlabel('Time [sec]')plt.grid()plt.legend()plt.subplots_adjust(hspace=0.35)#show the signal spetrum before and after filter.plt.figure(2);plt.subplot(2,1,1);yf=fft(data);xf=np.linspace(0.0,fs/2,n/2);plt.plot(xf,(np.abs(yf[0:n/2])));plt.subplot(2,1,2);yf_filter=fft(y);plt.plot(xf,(np.abs(yf_filter[0:n/2])));plt.grid();plt.show()


相關文章

聯繫我們

該頁面正文內容均來源於網絡整理,並不代表阿里雲官方的觀點,該頁面所提到的產品和服務也與阿里云無關,如果該頁面內容對您造成了困擾,歡迎寫郵件給我們,收到郵件我們將在5個工作日內處理。

如果您發現本社區中有涉嫌抄襲的內容,歡迎發送郵件至: info-contact@alibabacloud.com 進行舉報並提供相關證據,工作人員會在 5 個工作天內聯絡您,一經查實,本站將立刻刪除涉嫌侵權內容。

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.