功能描述：

      在信號處理中經常要把某些單頻(窄帶)干擾信號去除，例如求系統采集信號中的把工頻信號濾除。實際上有一個很好的方法，便是使用陷波器。在附件中給出了陷波器的設計技術，并舉了例子。

陷波器是無限沖擊響應(IIR)數字濾波器，該濾波器可以用以下常系數線性差分方程表示：

    由傳遞函數的零點和極點可以大致繪出頻率響應圖。在零點處，頻率響應出現極小值；在極點處，頻率響應出現極大值。因此可以根據所需頻率響應配置零點和極點，然后反向設計帶陷數字濾波器。考慮一種特殊情況，若零點在第1象限單位圓上，極點在單位圓內靠近零點的徑向上。為了防止濾波器系數出現復數，必須在z平面第4象限對稱位置配置相應的共軛零點、共軛極點。

clear all;
clc;


b=[1 -sqrt(2) 1];
a=[1 -sqrt(2)*0.999 0.999];
[db, mag, pha, grd, w]=freqz_m(b, a);        
subplot(221); plot(w*200/pi, db); title(' Magnitude Response' );
xlabel('frequency in Hz'); ylabel('dB'); axis([0, 100, -200, 5]);
set(gca, 'XTickMode', 'manual', 'XTick', [0, 50, 100]); 
set(gca, 'YTickmode', 'manual', 'YTick', [-200, -100, 0]); grid
title('Notch filter response');
t0=1:8000;
t=1:256;
t1=1:100;
t2=1:128;
x=sin(2*pi*50*t0/400)+0.5*sin(2*pi*100*t0/400);
x1=x(t);
y=filter(b,a,x1);
subplot(222); plot(x1);
title('Original waveform');
X=fft(x1);
subplot(223); plot(t2*400/256,abs(X(t2)));
xlabel('frequency in Hz'); ylabel('|H|'); axis([0, 200, 0, 150]);
title('Spectrum for original');
set(gca, 'XTickMode', 'manual', 'XTick', [0, 50, 100, 150]); 
set(gca, 'YTickmode', 'manual', 'YTick', [50, 100]); grid
y=filter(b,a,x);
x1=y(t+7600);
X=fft(x1);
subplot(224); plot(t2*400/256,abs(X(t2)));
xlabel('frequency in Hz'); ylabel('|H|'); axis([0, 200, 0, 150]);
title('Spectrum after filter'); 
set(gca, 'XTickMode', 'manual', 'XTick', [0, 50, 100, 150]); 
set(gca, 'YTickmode', 'manual', 'YTick', [50, 100]); grid


figure(2);
subplot(611);plot(x(t1)); axis([1, 100, -1.5, 1.5]); ylabel('input x');
set(gca, 'YTickmode', 'manual', 'YTick', [-1,-0.5,0, 0.5,1]); grid
subplot(612);plot(y); axis([1, 100, -1.5, 1.5]); ylabel('first');
set(gca, 'YTickmode', 'manual', 'YTick', [-1,-0.5,0,0.5,1]); grid
subplot(613);plot(y); axis([401, 500, -1.5, 1.5]); ylabel('second');
set(gca, 'YTickmode', 'manual', 'YTick', [-1,-0.5,0, 0.5,1]); grid
subplot(614);plot(y); axis([1201, 1300, -1.0, 1.0]); ylabel('forth');
set(gca, 'YTickmode', 'manual', 'YTick', [-0.5,0, 0.5]); grid
subplot(615);plot(y); axis([2000, 2100, -1.0, 1.0]); ylabel('sixth');
set(gca, 'YTickmode', 'manual', 'YTick', [-0.5,0, 0.5]); grid
subplot(616);plot(y); axis([3601, 3700, -1.0, 1.0]); ylabel('tenth');
set(gca, 'YTickmode', 'manual', 'YTick', [-0.5,0, 0.5]); grid

figure(3);
subplot(611);plot(y); axis([4401, 4500, -1, 1]); ylabel('twelfth');
set(gca, 'YTickmode', 'manual', 'YTick', [-0.5,0, 0.5]); grid
subplot(612);plot(y); axis([5201, 5300, -1.0, 1.0]); ylabel('fourteenth');
set(gca, 'YTickmode', 'manual', 'YTick', [-0.5,0, 0.5]); grid
subplot(613);plot(y); axis([6001, 6100, -1.0, 1.0]); ylabel('sixteenth');
set(gca, 'YTickmode', 'manual', 'YTick', [-0.5,0, 0.5]); grid
subplot(212);plot(y); axis([7601, 7650, -1.0, 1.0]); ylabel('twentieth');
set(gca, 'YTickmode', 'manual', 'YTick', [-0.5, 0, 0.5]); grid

function data=trapper(data,f,sf,n);
for i=1:n
fi=i*f;
omg=2*pi*fi/sf;
num=[1 -2*cos(omg) 1];
den=[1 -2*0.999*cos(omg) 0.999];
data=filter(num,den,data);
end
function [db, mag, pha, grd,w]=freqz_m(b,a);
%Modified version of freqz subroutine

[H,w]=freqz(b,a,1000,'whole');
H=(H(1:501))'; w=(w(1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(b,a,w);

end

fpga代做,fpga專業代做,simulink代做,simulink專業代做