matlab手动实现IIR滤波器

使用工具为matlab7.0
一，利用matlab fadtool工具获取SOS matrix 和 factor scales
1.在命令框中输入fdatool
2.配置滤波器信息，配置结果如下图

3.将SOS matrix 和 factor scales保存到workspace,保存的变量名为SOS和scale，并且生成my_filter.m文件，my_filter.m的文件内容如下
function Hd = my_filter
N = 3; % Order
Fc = 150; % Cutoff Frequency
h = fdesign.lowpass('N,Fc', N, Fc, Fs);
Hd = butter(h);

二.将SOS matrix 和 factor scales转化为A,B系数
1.在matlab命令对话框输入命令
g = prod(scale)
[b,a] = sos2tf(SOS,g)

三.手动编写代码实现IIR滤波器
IIR差分方程如下：

function F = filter()
% 采样频率
fs = 1000;
x = 0:1/fs:20;
% 构造的信号
sing = 0.5sin(2pi100x)+0.8cos(2pi180x);
% 将数据进行fft变换
function my_plot(data)
N = length(data);
xdft = fft(data);
xdft = xdft(1:N/2+1);
% psdx = abs(xdft).^2;
freq = 0:fs/N:fs/2;
plot(freq, (2/N)abs(xdft))
end
subplot(311)
my_plot(sing)
% A,B系数
b = [0.04953,0.14860,0.14860,0.04953];
a = [1,-1.16192,0.69594,-0.13776];
len = length(sing);
% 前几个时刻的输入值
pre_x = [0, 0, 0, 0];
% 前几个时刻的输出值
pre_y = [0, 0, 0, 0];
out = zeros(1,len);
for i=1:len
pre_x(1) = sing(i);
% 差分方程
out(i) = b(1)pre_x(1)+b(2)pre_x(2)+b(3)pre_x(3)+b(4)pre_x(4)-a(2)pre_y(2)-a(3)pre_y(3)-a(4)pre_y(4);
pre_y(1) = out(i);
for j = 4 : -1 : 2
pre_x(j) = pre_x(j-1);
pre_y(j) = pre_y(j-1);
end
end
subplot(312)
my_plot(out)
% 运行直接保存的.m文件
Hd = my_filter;
output = filter(Hd, sing);
subplot(313)
my_plot(output)
end

结果如下：