Matlab的：卷積帶通濾波器不切割不想要的頻率

我有在60kHz的採樣三個頻率分量的6ms的長的信號：Matlab的：卷積帶通濾波器不切割不想要的頻率

fs = 60000; 
T = 0.006; 
t = 0:1/fs:T; 
x = 0.3*sin(2*pi*2000*t) + sin(2*pi*5000*t) + 0.4*sin(2*pi*8000*t);

我與脈衝響應是兩者之間的差異的帶通濾波器正弦函數：

M = 151; 
N = 303; 
n = 0:(N-1); 
h = (sin(0.5760*pi*(n-M))-sin(0.3665*pi*(n-M)))./pi./(n-M); 
h(n==M) = 0.2094;

我設計與濾波器卷積的輸入的功能：

function y = fir_filter(h,x) 
y = zeros(1,length(x)+length(h)-1); 
for i = 1:length(x) 
for j = 1:length(h) 
    y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j); 
end 
end

然後應用的濾波器：

y = fir_filter(h,x);

這產生奇怪的結果：

figure(21) 
ax1 = subplot(311); 
plot(x); 
title('Input Signal'); 
ax2 = subplot(312); 
plot(h); 
title('FIR'); 
ax3 = subplot(313); 
plot(y); 
title('Output Signal'); 
linkaxes([ax1,ax2,ax3],'x') 
ax2.XLim = [0,length(y)];

由於過濾器是帶通，只有一個頻率分量預計生存。

我試着用yy = filter(h,1,[x,zeros(1,length(h)-1)]);和yyy = conv(h,x);得到了同樣的結果。請問任何人都可以解釋我做錯了什麼？謝謝！

來源

2015-09-15 brainkz

這是什麼情節？頻率響應？你的問題是fft鏡像？閱讀https://dsp.stackexchange.com/questions/4825/why-is-the-fft-mirrored或http://www.phys.nsu.ru/cherk/fft.pdf –

問題是，卷積我的FIR濾波器輸入產生錯誤的結果。結果應該接近純正弦曲線 – brainkz

您可以添加代碼來生成圖表嗎？ –

您的通帶不包括三個信號頻率分量中的任何一個。這可以直接在圖表上看到（第二個數字，包含脈衝響應，與信號相比變化太快）。或者，它可以從數字0.5760和0.3655在

h = (sin(0.5760*pi*(n-M))-sin(0.3665*pi*(n-M)))./pi./(n-M);

可以看出你爲什麼選擇這些號碼？信號的歸一化頻率爲[2 5 8]/60，即0.0333 0.0833 0.1333。它們都在通帶[.3665 .5760]以下。結果，濾波器極大地衰減了輸入信號的三個分量。假設你想隔離中心頻率分量（f = 5000Hz或f/fs = 0.08333標準化頻率）。你需要一個帶通濾波器，而不是讓頻率通過，並拒絕其他頻率。所以，你會使用，例如標準化的通帶[.06 .1]：

h = (sin(0.06*pi*(n-M))-sin(0.1*pi*(n-M)))./pi./(n-M); 
h(n==M) = (h(n==M+1)+h(n==M-1))/2; %// auto adjustment to avoid the 0/0 sample

與您的代碼的第二個問題是，它給出了兩個錯誤，因爲你指數h與0。爲了解決這個問題，改變

n = 0:(N-1);

到

n = 1:N;

和

y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j);

到

y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j+1);

因此，修改後的代碼，以隔離的中心頻率成分是：

fs = 60000; 
T = 0.006; 
t = 0:1/fs:T; 
x = 0.3*sin(2*pi*2000*t) + sin(2*pi*5000*t) + 0.4*sin(2*pi*8000*t); 

M = 151; 
N = 303; 
n = 1:N; %// modified 
h = (sin(0.06*pi*(n-M))-sin(0.1*pi*(n-M)))./pi./(n-M); %// modified 
h(n==M) = (h(n==M+1)+h(n==M-1))/2; %// modified 


y = zeros(1,length(x)+length(h)-1); 
for i = 1:length(x) 
for j = 1:length(h) 
    y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j+1); %// modified 
end 
end 

figure(21) 
ax1 = subplot(311); 
plot(x); 
title('Input Signal'); 
ax2 = subplot(312); 
plot(h); 
title('FIR'); 
ax3 = subplot(313); 
plot(y); 
title('Output Signal'); 
linkaxes([ax1,ax2,ax3],'x') 
ax2.XLim = [0,length(y)];

結果如下。