-
本節會繪製七個圖
第一個子圖為100Hz的正弦波
第二列的兩個子圖為附加不同鎮雜訊的100Hz的正弦波
第三列的兩個子圖為第二列的兩個子圖的Welch功率譜
第四列的兩個子圖為第二列的兩個子圖的傅立葉頻域譜
# -*- coding: UTF-8 -*- import numpy as np import scipy as sp import scipy.signal as signal import matplotlib.pyplot as plt from scipy.fftpack import fft fs = 1e4 #sampleing rate N = 5e2 #sampleing number amp = 1 #ampitude freq = 100.0 #frequency noise_power = 0.001 * fs / 2 time = np.arange(N) / fs x = amp * np.sin(2*np.pi*freq*time) #x2 = x + np.random.normal(amp*0.1, size=time.shape) noise1 = sp.randn(time.size) noise0p5 = 0.5*sp.randn(time.size) x1 = x + noise1 x2 = x + noise0p5 #Estimate power spectral density using Welch’s method. #Welch’s method computes an estimate of the power spectral density by dividing the data into overlapping segments, computing a modified #periodogram for each segment and averaging the periodograms. fff = np.zeros(6) fig, (fff) = plt.subplots(4,2) f, Pwelch_spec = signal.welch(x1, fs, scaling='spectrum') f, Pwelch_spec0p5 = signal.welch(x2, fs, scaling='spectrum') y = sp.fftpack.fft(x1,N)/N y0p5 = sp.fftpack.fft(x2,N)/N print(y.size) f2 = (fs/N)*sp.arange(N/2+1) print(f2.size) fff[0,0].plot(time, x) fff[1,0].plot(time, x1) fff[2,0].plot(f, Pwelch_spec) fff[3,0].plot(f, Pwelch_spec0p5) fff[1,1].plot(time, x2) fff[2,1].plot(f2, 2*abs( y[:int(N/2)+1] )) fff[3,1].plot(f2, 2*abs( y0p5[:int(N/2)+1] )) plt.show()
# -*- coding: UTF-8 -*- import numpy as np import scipy as sp import scipy.signal as signal import matplotlib.pyplot as plt from scipy.fftpack import fft fs = 1e4 #取樣率 N = 5e2 #取樣數 amp = 1 #振幅 freq = 100.0 #頻率 noise_power = 0.001 * fs / 2 ＃雜訊強度 time = np.arange(N) / fs ＃每個取樣的時間 x = amp * np.sin(2*np.pi*freq*time) ＃產生正弦波 noise1 = sp.randn(time.size) ＃利用randn給出N個依照normal distribution 的亂數 noise0p5 = 0.5*sp.randn(time.size) sp.randn(time.size) ＃利用randn給出N個依照normal distribution 的亂數, 但最大振幅只有0.5 x1 = x + noise1 ＃產生一個帶雜訊的正弦波(雜訊最大振幅1) x2 = x + noise0p5 ＃產生一個帶雜訊的正弦波(雜訊最大振幅0.5) #Estimate power spectral density using Welch’s method. #Welch’s method computes an estimate of the power spectral density by dividing the data into overlapping segments, computing a modified #periodogram for each segment and averaging the periodograms. #fff = np.zeros(6) #利用zero函數產生六個0的數列, 作為圖形的引數 fig, fff = plt.subplots(4,2) f, Pwelch_spec = signal.welch(x1, fs, scaling='spectrum') f, Pwelch_spec0p5 = signal.welch(x2, fs, scaling='spectrum') y = sp.fftpack.fft(x1,N)/N y0p5 = sp.fftpack.fft(x2,N)/N print(y.size) f2 = (fs/N)*sp.arange(N/2+1) print(f2.size) fff[0,0].plot(time, x) fff[1,0].plot(time, x1) fff[2,0].plot(f, Pwelch_spec) fff[3,0].plot(f, Pwelch_spec) fff[3,0].plot(f, Pwelch_spec0p5) fff[1,1].plot(time, x2) fff[2,1].plot(f2, 2*abs( y[:N/2+1] )) fff[3,1].plot(f2, 2*abs( y[:N/2+1] )) fff[3,1].plot(f2, 2*abs( y0p5[:N/2+1] )) plt.show()
