我认为，现在滤波信号的 FFT 图上噪声区域的幅度应该更低。可能是scipy.fftpack.lfilter().
我试图在现有信号中添加噪声，但没有结果。
为什么噪声区域中滤波信号（绿色）的 FFT 幅度如此之高？
更新：
FFT 幅度的 300 dB 是非物理的 - 很明显，这是由于 Python 环境中的 64 位浮点数造成的。
滤波信号（绿色）的 FFT 具有如此低的 dB（约 67 dB），因为所有信号都有约 4000 个样本（不是图片上的“每秒”，小错误，不重要），采样率 = 200 个样本/秒。1 个频率分档=200/4000/2=0.025Hz，显示 2000 个分档。
如果我们采用更长的信号，我们将获得更高的每个 bin 的频率分辨率（即 40 000 个样本，1 freq bin = 200/40 000/2 = 0.0025 Hz）。我们还得到了滤波信号的 FFT ~87 dB。
（数字 67 dB 和 87 dB 似乎是非物理的，因为初始信号 SNR 为 300dB，但我尝试在现有信号中添加一些噪声并获得相同的数字）
如果您想获得与信号中样本数无关的 FFT 图像，则应使用加窗 FFT 和滑动窗口来计算全信号 FFT。
'''
Created on 13.02.2013, python 2.7.3
@author:
'''
from numpy.random import normal
from numpy import sin, pi, absolute, arange, round
#from numpy.fft import fft
from scipy.fftpack import fft, ifft
from scipy.signal import kaiserord, firwin, lfilter, freqz
from pylab import figure, clf, plot, xlabel, ylabel, xlim, ylim, title, grid, axis, show, log10,\
subplots_adjust, subplot
def filter_apply(filename):
pass
def sin_generator(freq_hz = 1000, sample_rate = 8000, amplitude = 1.0, time_s = 1):
nsamples = round(sample_rate * time_s)
t = arange(nsamples) / float(sample_rate.__float__())
signal = amplitude * sin(2*pi*freq_hz.__float__()*t)
return signal, nsamples
def do_fir(signal, sample_rate):
return signal
#-----------------make a signal---------------
freq_hz = 10.0
sample_rate = 400
amplitude = 1.0
time_s = 10
a1, nsamples = sin_generator(freq_hz, sample_rate, amplitude, time_s)
a2, nsamples = sin_generator(50.0, sample_rate, 0.5*amplitude, time_s)
a3, nsamples = sin_generator(150.0, sample_rate, 0.5*amplitude, time_s)
mu, sigma = 0, 0.1 # mean and standard deviation
noise = normal(mu, sigma, nsamples)
signal = a1 + a2 + a3 # + noise
#----------------create low-pass FIR----
# The Nyquist rate of the signal.
nyq_rate = sample_rate / 2.0
# The desired width of the transition from pass to stop,
# relative to the Nyquist rate. We'll design the filter
# with a 5 Hz transition width.
width = 5.0/nyq_rate
# The desired attenuation in the stop band, in dB.
ripple_db = 60.0
# Compute the order and Kaiser parameter for the FIR filter.
N, beta = kaiserord(ripple_db, width)
print 'N = ',N, 'beta = kaiser param = ', beta
# The cutoff frequency of the filter.
cutoff_hz = 30.0
# Use firwin with a Kaiser window to create a lowpass FIR filter.
# Length of the filter (number of coefficients, i.e. the filter order + 1)
taps = firwin(N, cutoff_hz/nyq_rate, window=('kaiser', beta))
# Use lfilter to filter x with the FIR filter.
filtered_signal = lfilter(taps, 1.0, signal)
#----------------plot signal----------------------
hh,ww=2,2
figure(figsize=(12,9))
subplots_adjust(hspace=.5)
#figure(1)
subplot(hh,ww,1)
# existing signal
x = arange(nsamples) / float(sample_rate)
# The phase delay of the filtered signal.
delay = 0.5 * (N-1) / sample_rate
# original signal
plot(x, signal, '-bo' , linewidth=2)
# filtered signal shifted to compensate for
# the phase delay.
plot(x-delay, filtered_signal, 'r-' , linewidth=1)
# Plot just the "good" part of the filtered signal.
# The first N-1 samples are "corrupted" by the
# initial conditions.
plot(x[N-1:]-delay, filtered_signal[N-1:], 'g', linewidth=2)
xlabel('time (s)')
ylabel('amplitude')
axis([0, 1.0/freq_hz*2, -(amplitude*1.5),amplitude*1.5]) # two periods of freq_hz
title('Signal (%d samples)' % nsamples)
grid(True)
#-------------- FFT of the signal
subplot(hh,ww,2)
signal_fft=fft(signal)
filtered_fft =fft(filtered_signal[N-1:])
# existing signal
y = 20*log10( ( abs( signal_fft/nsamples )*2.0)/max( abs( signal_fft/nsamples )*2.0) )# dB Amplitude
x = arange(nsamples)/float(nsamples)*float(sample_rate)
# filtered signal
y_filtered = 20*log10( (abs(filtered_fft/ (nsamples - N + 1) )*2.0)/max(abs(signal_fft/ (nsamples - N + 1) )*2.0) )# dB Amplitude
x_filtered = arange(nsamples - N + 1)/float(nsamples - N + 1)*float(sample_rate)
yy = fft(ifft(filtered_fft))
plot(x,y, linewidth=1)
plot(x_filtered, y_filtered, 'g', linewidth=2)
xlim(0, sample_rate/2) # compensation of mirror (FFT imaginary part)
xlabel('freq (Hz)')
ylabel('amplitude, (dB)')
title('Signal (%d samples)' % nsamples)
grid(True)
#--------------FIR ampitude response
subplot(hh,ww,3)
w, h = freqz(taps, worN=8000)
#plot((w/pi)*nyq_rate, absolute(h), linewidth=2)
plot((w/pi)*nyq_rate, 20*log10(absolute(h)/1.0),'r', linewidth=1)
xlabel('Frequency (Hz)')
#ylabel('Gain -blue')
ylabel('Gain (dB)')
title('Frequency Response')
#ylim(-0.05, 1.05)
grid(True)
#--------------FIR coeffs
subplot(hh,ww,4)
plot(taps, 'bo-', linewidth=2)
title('Filter Coefficients (%d taps)' % N)
grid(True)
show()