现有（合成）信号与滤波信号之间的 FFT 幅度差异很大

Question

我认为，过滤后信号的 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()

Answer 1

我认为你 -300dB 的噪音相当低。请记住，这是一个对数刻度，所以我们说的是 64 位分辨率下的几个数字。

使用双精度浮点数（64 位）你不会得到任何更低的值。