0

フィルター処理された信号の FFT プロット上のノイズ領域の振幅は、今よりも低くなければならないと思います。おそらく、小さな数値偏差/誤差がscipy.fftpack.lfilter().

既存の信号にノイズを追加しようとしましたが、結果はありません。

ノイズ領域でフィルタリングされた信号 (緑) の FFT 振幅が非常に高いのはなぜですか?

ここに画像の説明を入力

更新:
FFT 振幅の 300 dB は非物理的です。これは明らかです。Python 環境の 64 ビット浮動小数点が原因です。

フィルター処理された信号 (緑) の FFT は、すべての信号で ~4000 サンプル (写真のように「1 秒あたり」ではない、小さな間違い、重大ではない) のため、非常に低い dB (約 67 dB) であり、サンプルレート = 200 サンプル/秒 1 周波数ビン = 200/4000/2=0.025Hz、2000 ビンが表示されます。

より長い信号を取得すると、ビンごとにより多くの周波数分解能が得られます (つまり、40 000 サンプル、1 周波数ビン = 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()
4

1 に答える 1

2

-300dBのノイズはかなり低いと思います。これは対数スケールであるため、64 ビットの解像度で数桁について話していることに注意してください。

倍精度浮動小数点数 (64 ビット) を使用すると、それ以上低くなりません。

于 2013-02-18T12:52:13.823 に答える