问题描述
我认为，滤波信号的FFT曲线上的噪声区域中的幅度应该更低，然后现在. scipy.fftpack.lfilter()中可能出现较小的数字偏差/误差.
I think, what amplitude in the noised regions on the FFT-plot of the filtered signal should be lower, then now. Probably, small numeric deviations/errors presented in the scipy.fftpack.lfilter().
我试图将噪声添加到现有信号中，但没有结果.
I tried to add noise to the existing signal, but no result.
为什么噪声区域中的滤波信号(绿色)的FFT幅度这么高?
Why FFT-amplitude of the filtered signal (green) in the noise regions are so high?
更新:
300 dB的FFT幅度是非物理的-显而易见，这是由于Python环境中的64位浮点数引起的.
Update:
300 dB of an FFT-amplitude are non-physical - it's clear, it's due to 64bit float in the Python-environment.
由于在所有信号上进行了约4000个采样(并非图片中的每秒"，几乎没有错误，不是很关键)，所以滤波信号(绿色)的FFT具有如此低的dB(约为67 dB)，采样率= 200个样本/秒1个频率bin = 200/4000/2 = 0.025Hz，显示了2000个bin.
FFT of the filtered signal (green) has so low dB (as ~67 dB) due to ~4000 samples at all signal (not 'per second' as on the picture, little mistake, not critical), sample rate = 200 samples/sec. 1 frequency bin=200/4000/2=0.025Hz, 2000 bins shown.
如果采用更长的信号，则每个bin的频率分辨率更高(即40 000个样本，1 freq bin = 200/40 000/2 = 0.0025 Hz).而且我们还获得了滤波后的信号〜87 dB的FFT.
If we take a longer signal, we get a more frequency resolution per bin (i.e. 40 000samples, 1 freq bin = 200/40 000/2 = 0.0025 Hz). And also we get FFT of a filtered signal ~87 dB.
(数字67 dB和87 dB似乎是非物理的，因为初始信号的SNR为300dB，但我尝试在现有信号上添加一些噪声并获得相同的数字)
(Numbers 67 dB and 87 dB seems to be non-physical, because of an initial signal SNR 300dB, but I try to add some noise to the existing signal and get the same numbers)
如果要获取与信号中的样本数量无关的FFT图片，则应使用加窗FFT和滑动窗口来计算全信号FFT.
If you want to get a picture of a non-depend FFT to the number of samples in a signal, you should use a windowed FFT and slide window to compute a whole-signal-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()
推荐答案
我认为-300dB的噪声非常低.请记住，这是一个对数刻度，因此我们在谈论64位分辨率下的几个数字.
I think you -300dB of noise is pretty low. Keep in mind the this is a log scale, so we are talking about a couple of digits at 64bit resolution.
使用双精度浮点数(64位)，您将获得更低的精度.
Using double precision float (64bit) you won't get any lower.
这篇关于现有(合成)信号与滤波后信号之间的FFT幅度差异较大的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！