设计巴特沃斯滤波器实现高通滤波

首先在MATLAB上通过设计一个巴特沃斯滤波器来实现一个高通滤波器。

%设计一个巴特沃斯高通滤波器
%滤除一个300HZ的正弦信号
clear;
fn = 10000;%采样频率
fp = 900;%通带截止频率
fs = 600;%阻带截止频率
Rp = 2;%通带最大衰减
Rs = 20;%阻带最小衰减
figure(1);
Wp = fp/(fn/2);
Ws = fs/(fn/2);
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[b,a] = butter(n,Wn,'high');

[H,F] = freqz(b,a,900,10000);%数字滤波器的频率响应
subplot(211);
plot(F,20*log10(abs(H)))
axis([0 4000 -30 3]);

subplot(212);
pha = angle(H)*180/pi;
plot(F,pha);
% axis([0 8000 -200 200]);
grid on;
figure(2);
f1 = 2000;
f2 = 300;
t = 0:1/fn:1;
x = sin(2*pi*f1*t)+cos(2*pi*f2*t)/6;
y = filter(b,a,x);
subplot(211);
Y = fft(x);
n=0:length(x)-1;
plot(n,Y);
subplot(212);
Y1 = fft(y);
n=0:length(y)-1;
plot(n,Y1);
grid on;
figure(4);
plot(t,x);
subplot(111);
figure(5);
plot(t,y);
subplot(111);

结果如下：
滤波器的频率特性：

在600HZ时-20分贝
输入波形：

可以从幅度看出输入波形为sin+cos/6；
输出波形：

输入和输出的频谱：

观察滤波器【b,a】[b,a] = butter(n,Wn,‘high’);
a = 1 -4.67625140576819 9.64877723309598 -11.3207795685212 8.12672213051155 -3.55976931154069 0.879191790281082 -0.0942959539620253
b = 0.307076464013131 -2.14953524809191 6.44860574427574 -10.7476762404596 10.7476762404596 -6.44860574427574 2.14953524809191 -0.307076464013131
则系统函数 y(n)= 4.67625140576819y(n-1)-9.64877723309598y(n-2)…+0.0942959539620253y(n-7)+0.307076464013131x(n)-2.14953524809191x(n-1)…-0.307076464013131x(n-7)
所以通过这个编写CCS的代码：

#include"math.h"

#define IIRNUMBER 8
#define SIGNAL1F 1000
#define SIGNAL2F 4500
#define SAMPLEF  10000
#define PI 3.1415926

float InputWave();
float IIR();

float fBn[IIRNUMBER]={0,4.67625140576819,-9.64877723309598,11.3207795685212,-8.12672213051155,3.55976931154069,-0.879191790281082,0.0942959539620253};
float fAn[IIRNUMBER]={0.307076464013131,-2.14953524809191,6.44860574427574,-10.7476762404596,10.7476762404596,-6.44860574427574,2.14953524809191,-0.307076464013131};
float fXn[IIRNUMBER]={ 0.0 };
float fYn[IIRNUMBER]={ 0.0 };
float fInput,fOutput;
float fSignal1,fSignal2;
float fStepSignal1,fStepSignal2;
float f2PI;
int i;
float fIn[256],fOut[256];
int nIn,nOut;

main()
{
	nIn=0; nOut=0;
	fInput=fOutput=0;
	f2PI=2*PI;
	fSignal1=0.0;
	fSignal2=PI*0.1;//初始相位
	fStepSignal1=2*PI*300;
	fStepSignal2=2*PI*2000;
	while ( 1 )
	{
		fInput=InputWave();
		fIn[nIn]=fInput;
		nIn++; nIn%=256;
		fOutput=IIR();
		fYn[0]=fOutput;
		fOut[nOut]=fOutput;
		nOut++;				// break point
		if ( nOut>=256 )
		{
			nOut=0;		
		}
	}
}
/*生成采样率10000HZ，频率300HZ的正弦信号+1/6频率2000HZ的余弦信号*/
float InputWave()
{
	for ( i=IIRNUMBER-1;i>0;i-- )
	{
		fXn[i]=fXn[i-1];
		fYn[i]=fYn[i-1];
	}
	fXn[0]=sin((double)fSignal1)+cos((double)fSignal2)/6.0;

	fSignal1+=fStepSignal1/SAMPLEF;
	if ( fSignal1>=f2PI )	fSignal1-=f2PI;
	fSignal2+=fStepSignal2/SAMPLEF;
	if ( fSignal2>=f2PI )	fSignal2-=f2PI;
	return(fXn[0]);
}
//求y(n)=4.676251y(n-1)+....++0.0942959y(n-7)+0.3070764x(n)+....-0.307076x(n-7)
float IIR()
{
	float fSum;
	fSum=0.0;
	for ( i=0;i<IIRNUMBER;i++ )
	{
		fSum+=(fXn[i]*fAn[i]);
		fSum+=(fYn[i]*fBn[i]);

	}
	return(fSum);
}

debug结果如下：
输入信号：

滤波后输出信号：

输入信号频谱：

输出信号频谱：