现有信号与滤波信号的FFT幅度差异大

Question

我觉得，在经过滤波的信号的FFT图上，噪声区域的幅度应该比现在低。可能是因为在使用scipy.fftpack.lfilter()时出现了一些小的数值偏差或错误。

我尝试在现有信号中添加噪声，但没有得到任何结果。

为什么经过滤波的信号（绿色）的FFT幅度在噪声区域会这么高呢？

在这里输入图片描述

更新：
300 dB的FFT幅度是非物理的，这很明显，是因为在Python环境中使用了64位浮点数。

经过滤波的信号（绿色）的FFT幅度很低（大约67 dB），这是因为整个信号只有大约4000个样本（不是像图中所示的“每秒”那样，稍微有点错误，但不重要），采样率为200样本/秒。一个频率区间=200/4000/2=0.025Hz，显示了2000个区间。

如果我们使用更长的信号，就能在每个区间获得更好的频率分辨率（比如说40000个样本，1个频率区间=200/40000/2=0.0025 Hz）。这样，我们得到的经过滤波的信号的FFT幅度大约是87 dB。

（67 dB和87 dB的数字看起来是非物理的，因为初始信号的信噪比是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()

信号处理采样率数值稳定性窗口函数 fft 滤波器噪声频率分辨率

现有信号与滤波信号的FFT幅度差异大

1 个回答

撰写回答