matlab语言与应用 08 数值插值、函数逼近

现代科学运算-matlab语言与应用

$\qquad \qquad$ 东北大学 http://www.icourses.cn/home/ (薛定宇)
《高等应用数学问题的MATLAB求解（第四版）》
代码可在matlab r2016b 运行。

08 数据插值、函数逼近问题的计算机求解

08.01 一维插值问题

一维插值 interp1() 函数的调用格式为
$y_1 = \mathbb{interp1}(x, y, x_1, method)$
插值方法：
Linear – 默认方法
Pchip – 三次函数插值
Nearest – 最近点插值
Spline – 样条插值，推荐使用

运用外推法，再区间[$x_1, x_N$] 外的点的插值
$y_1 = \mathbb{interp1}(x, y, x_1, method, 'extrap')$

例8-1 数据插值
已知的数据点来自函数
$f(x) = (x^2 - 3x + 5)e^{-5x} \sin x$
生成的数据进行插值处理
得出较平滑的插值曲线
生成的数据显示

x = 0: .12: 1;
y = (x.^2 - 3*x + 5).*exp(-5*x).*sin(x);
plot(x, y, x, y, 'o');

不同的插值算法

x1 = 0: .02: 1;
y1 = interp1(x, y, x1);
y0 = (x1.^2 - 3*x1 + 5).*exp(-5*x1).*sin(x1);
y2 = interp1(x, y, x1, 'pchip');
y3 = interp1(x, y, x1, 'spline');
y4 = interp1(x, y, x1, 'nearest');
plot(x1, [y1', y2', y3', y4'], ':', x, y, 'o', x1, y0)
% 验证结果
[max(abs(y0(1:49) - y2(1:49))), max(abs(y0 - y3)), max(abs(y0 - y4))]

例8-2 徒手绘制规划曲线
用户可以在 x-y 平面内选择一些点
x轴排序，密集点插值，绘图

function sketcher(vis)
  x = [];
  y = [];
  i = 1;
  h = [];
  axes(gca);
  while 1
    [x0,y0,but] = ginput(1); 
    if but == 1
      x = [x, x0];
      y = [y, y0]; 
      h(i) = line(x0, y0);
      set(h(i), 'Marker', 'o');
      i = i + 1;
    else
      break 
    end
  end
  if nargin == 1
    delete(h);
  end 
  [x ii] = sort(x);
  y = y(ii);
  xx = [x(1):(x(end)-x(1))/100: x(end)]; 
  yy = interp1(x, y, xx, 'spline');
  line(xx, yy) 

Lagrange插值算法及应用
$已知x_i, y_i 点, i = 1, 2, \cdots, N$
生成一个插值x向量
插值算法为 $\phi(x) = \sum \limits_{i=1}^{N} y_i \prod \limits_{j = 1, j \neq i}^{N} \dfrac{x - x_j}{(x_i - x_j)}$
matlab函数调用格式： y = lagrange($x_0, y_0, x$)

例8-3 Runge现象
测试函数 $f(x) = 1/(1+25x^2), -1 \leq x \leq 1$
生成数据，进行Lagrange插值

function y = lagrange(x0, y0, x)
  ii = 1:length(x0);
  y=zeros(size(x));
  for i = ii
    ij = find(ii ~= i);
    y1 = 1;
    for j = 1:length(ij)
      y1 = y1.*(x-x0(ij(j)));
    end
    y = y+y1*y0(i)/prod(x0(i)-x0(ij));
  end

x0 = -1+2*[0:10]/10;
y0 = 1./(1+25*x0.^2);
x = -1: .01:1;
y = lagrange(x0, y0, x);
ya = 1./ (1+25*x.^2);
plot(x, ya, x, y, ':');

调用interp1()函数

y1 = interp1(x0, y0, x, 'pchip');
y2 = interp1(x0, y0, x, 'spline');
plot(x, ya, x, y1, ':', x, y2, '--');

例8-4 人口预测
某省人口数据在census.xls给出
可以根据给出的数据预报未来的人口

X = xlsread('E:\\Work\\201806\\math_4th\\census.xls', 'B5:C67');
t = X(:, 1);
p = X(:, 2);
t0 = t(1:5:end);
p0 = p(1:5:end);
t1 = 1949:2015;
y = interp1(t0, p0, t1, 'spline', 'extrap');
plot(t, p, t1, y, t0, p0, 'o');

已知样本点的定积分计算
已知函数散点$x_0, y_0$

function y = quadspln(x0, y0, a, b)
  f = @(x) interp1(x0, y0, x, 'spline');
  y = integral(f, a, b);

函数调用格式: I = quadspln(x0, y0, a, b)

例8-5 利用插值求数值积分
利用样条插值三分求解 $\int _0^{\pi} \sin x d x$
用30个采样点求解:

x0 = 0: pi/30: pi;
y0 = sin(x0);
I = quadspln(x0, y0, 0, pi)

08.02 二维与高维插值问题

二维网格数据的插值问题
样本点是网格形式给出的
二维插值的函数
$z_1 = interp2(x_0, y_0, z_0, x_1, y_1, method)$
其中，
$x_0, y_0, z_0 为已知二维网格数据x_1, y_1$
$x, y 为要插值点的网格矩阵z_1$
选项类似于 interp1，建议使用 spline 选项
这里的数据一定要采用网格数据
网格数据应该由 meshgrid函数生成

例8-7 网格数据的插值举例
根据下述函数生成一些稀疏的网格数据
由下面的函数可以生成样本点网格数据
$z = f(x, y) = (x^2 - 2x)e^{-x^2-y^2-xy}$
进行各种插值拟合，并比较拟合结果
绘制已知数据的网格图:

[x, y] = meshgrid(-3: .6 : 3, -2: .4 : 2);
z = (x.^2 - 2*x) .* exp(-x.^2-y.^2-x.*y);
surf(x, y, z);
axis([-3, 3, -2, 2, -0.7, 1.5]);

不同插值算法的比较
默认插值算法进行插值：

[x1, y1] = meshgrid(-3: .2: 3, -2: .2: 2);
z1 = interp2(x, y, z, x1, y1);
surf(x1, y1, z1);
axis([-3, 3, -2, 2, -0.7, 1.5]);

三次插值和样条插值，并分析误差

[x, y] = meshgrid(-3: .6 : 3, -2: .4 : 2);
z = (x.^2 - 2*x) .* exp(-x.^2-y.^2-x.*y);
[x1, y1] = meshgrid(-3: .2: 3, -2: .2: 2);
z1 = interp2(x, y, z, x1, y1, 'cubic');
z2 = interp2(x, y, z, x1, y1, 'spline');
surf(x1, y1, z1);
figure;
surf(x1, y1, z2);
z = (x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
figure;
surf(x1, y1, abs(z-z1));
figure;
surf(x1, y1, abs(z-z2));

二维一般分布数据的插值问题
实际应用中能测到的并不总是网格分布的数据，经常需要解决样本点为散点的插值
griddata()函数的调用格式：z = griddata($x_0, y_0, z_0, x, y$, ‘v4’)
其中, $x_0, y_0, z_0$是已知的二维样本点向量
x,y 是期望的插值位置
z表示插值的结果，维数和x,y一致
‘v4’是matlab 4.0 版本中提供的插值算法

例8-8 随机分布样本点的插值
给定函数 $z = f(x, y) = (x^2 - 2x)e^{-x^2-y^2-xy}$
在矩形区域$x \in [-3, 3], y \in [-2, 2]$内随机生成一组样本点$z_i$
用griddata()进行插值处理, 并误差分析
数据生成与显示

x = -3 + 6*rand(200, 1);
y = -2 + 4*rand(200, 1);
z = (x.^2 - 2*x).*exp(-x.^2-y.^2-x.*y);
plot(x, y, 'x');
figure;
plot3(x, y, z, 'x');

不同插值方法比较
使用’cubic’和’v4’算法：

[x1, y1] = meshgrid(-3: .2: 3, -2: .2: 2);
z1 = griddata(x, y, z, x1, y1, 'cubic');
surf(x1, y1, z1);
z2 = griddata(x, y, z, x1, y1, 'v4');
figure;
surf(x1, y1, z2);
% 比较
figure;
z0 = (x1.^2 - 2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
surf(x1, y1, abs(z0-z1));
figure;
surf(x1, y1, abs(z0-z2));

例8-9 剔除某些样本点
给定的样本点在x-y平面分布较均匀，现在人为剔除某些点,表明已知数据分布不均匀，这时再进行插值分析,观察插值效果。
剔除圆内样本点

x = -3 + 6*rand(200, 1);
y = -2 + 4*rand(200, 1);
ii = find((x+1).^2 + (y+0.5).^2 > 0.5^2);
x = x(ii);
y = y(ii);
z = z(ii);
plot(x, y, 'x');
t = [0: .1: 2*pi, 2*pi];
x0 = -1 + 0.5*cos(t);
y0 = -0.5 + 0.5*sin(t);
line(x0, y0);

重新插值

% 用新样本点拟合曲面
[x1, y1] = meshgrid(-3: .2: 3, -2: .2: 2);
z1 = griddata(x, y, z, z1, y1, 'v4');
surf(x1, y1, z1);
% 误差分析
z0 = (x1.^2 - 2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
surf(x1, y1, abs(z0-z1)),
axis([-3, 3, -2, 2, 0, 0.1]);
% 误差等高线图
figure;
contour(x1, y1, abs(z0-z1), 30);
hold on;
plot(x, y, 'x');
hold off;
line(x0, y0)

高维插值问题
三维的网格数据生成：
[x, y, z] = meshgrid($x_1, y_1, z_1$);
高维网格数据的生成
$[x_1, \cdots, x_n] = \mathbb{ndgrid}(v_1, \cdots, v_n);$
插值函数
网格数据插值 interpn
散点数据插值 griddatan

数据格式转换
ndgrid与meshgrid格式不同

% 格式转换方法
function [x1, y1, z1, v1] = mesh2nd(x, y, z, v)
  switch nargin
    case 3, x1 = x.'; y1 = y.'; z1 = z.';
    case 4, z1 = z;
      for i = 1: size(x, 3), x1(:, :, i) = x(:, :, i).';
        y1(:, :, i) = y(:, :, i).';
        v1(:, :, i) = v(:, :, i).';
      end;
    otherwise error('Error in input arguments')
  end;

例8-10 三元函数插值
假设已知某三元函数
$V(x, y, z) = \sqrt{x^x + y^{(x+y)/2} + z^{(x+y+z)/3}}$
可以通过该函数生成一些网格样本点，是根据样本点进行拟合，并给出拟合误差

[x, y, z] = meshgrid(0: 0.3: 2);
[x0, y0, z0] = meshgrid(0: 0.1: 2);
V = sqrt(x.^x + y.^((x+y)/2) + z.^((x+y+z)/3));
V0 = sqrt(x0.^x0 + y0.^((x0+y0)/2) + z0.^((x0+y0+z0)/3));
V1 = interp3(x, y, z, V, x0, y0, z0, 'spline');
vol_visual4d(x0, y0, z0, V1);

08.03 样条插值问题

三次样条函数及其matlab表示
样本点$(x_i, y_i) \; (i = 1, 2, \cdots, n)$
其中$x_1 < x_2 < \cdots < x_n$
S(x) 为三次样条函数的三个条件：
$S(x_i) = y_i (i = 1, 2, \cdots, n)$
每个子区间$[x_i, x_{i+1}]$ 上, S(x)位三次多项式
$S(x) = c_{i1}(x - x_i)^3 + c_{i2}(x - x_i)^2 + c_{i3}(x - x_i) + c_{i4}$
S(x) 在整个区间 [x_i, x_{i+1}] 上有连续1-2阶导数

例8.12 正弦函数插值
稀疏数据的三次样条插值结果
$I = \int_0^{\pi} \sin t dt$

x0 = [0, 0.4, 1 2, pi];
y0 = sin(x0);
sp = csapi(x0, y0);
fnplt(sp, ':');
hold on;
ezplot('sin(t)', [0, pi]);
plot(x0, y0, 'o');

多个自变量网格数据插值
网格数据三次样条插值类:
S = csapi($\lbrace x_1, x_2, \cdots, x_n \rbrace$, z)
$x_i$为自变量的网格数据
$z_i$为网格数据的样本点
得出的S是三次样条函数对象
得出的样条函数类可用下面函数处理
fnplt(), fnval()

例8-14 二维数据插值
给定$z = f(x, y) = (x^2-2x)e^{-x^2-y^2-xy}$
用三次样条插值方法得出网格数据的样条插值拟合，并绘制曲面

x0 = -3: .6: 3;
y0 = -2: .4: 2;
[x, y] = ndgrid(x0, y0);
z = (x.^2 - 2*x) .* exp(-x.^2 - y.^2 -x.*y);
sp = csapi({x0, y0}, z);
fnplt(sp);

B样条函数及其matlab表示
建立B样条函数对象
S = spapi(k, z, y)
其中k为用户选定的B样条阶次，一般选择k= 4, 5

例8-15 插值类的比较
给定正弦函数 $f(x) = \sin x$
生成稀疏数据
进行5次B样条插值函数拟合
与三次样条函数拟合的结果比较

x0 = [0, 0.4, 1 2, pi];
y0 = sin(x0);
ezplot('sin(t)', [0, pi]);
hold on;
sp1 = csapi(x0, y0);
fnplt(sp1, '--');
sp2 = spapi(5, x0, y0);
fnplt(sp2, ':');
hold off;

另一个函数
原函数 $f(x) = (x^2-3x+5)e^{-5x} \sin x$
插值效果比较

x = 0: .12: 1;
y = (x.^2-3*x+5).*exp(-5*x).*sin(x);
ezplot('(x^2-3*x)*exp(-5*x)*sin(x)', [0, 1]);
hold on;
sp1 = csapi(x, y);
fnplt(sp1, '--');
sp2 = spapi(5, x, y);
fnplt(sp2, ':');
hold off;

08.04 基于样条插值的数值微积分

基于样条插值的数值微分运算
求单变量k阶导数的函数调用格式
$S_d$ = fnder(S, k)
多变量函数的偏导数的函数调用格式
$S_d$ = fnder(S, [$k_1, \cdots, k_n$])
S可用三次或B样条类

例8-16数值微分
给定$f(x) = (x^2-3x+5)e^{-5x} \sin x$
生成一些数据点
用三次分段多项式样条函数与B样条插值函数，求出该函数的导数

syms x;
f = (x^2-3*x+5)*exp(-5*x)*sin(x);
ezplot(diff(f), [0, 1]);
hold on;
x = 0: .12: 1;
y = (x.^2-3*x+5).*exp(-5*x).*sin(x);
% 数值微分
figure;
sp1 = csapi(x, y);
dsp1 = fnder(sp1, 1);
fnplt(dsp1, '--');
sp2 = spapi(5, x, y);
dsp2 = fnder(sp2, 1);
fnplt(dsp2, ':');
hold off;
% 高阶导数，如三阶导数
figure;
ezplot(diff(f, 3), [0, 1]);
hold on;
dsp1 = fnder(sp1, 3);
fnplt(dsp1, '--');
dsp2 = fnder(sp2, 3);
fnplt(dsp2, ':');
hold off;

例8-17 二元函数
给定 $z = f(x, y) = (x^2 - 2x)e^{-x^2-y^2-xy}$
生成一些数据点
利用数值插值的方法拟合曲面 $\dfrac{partial ^z }{\partial x \partial y}$
并与解析解法绘制出的曲面相比较

x0 = -3: .3 : 3;
y0 = -2: .2: 2;
[x, y] = ndgrid(x0, y0);
z = (x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
sp = spapi({5, 5}, {x0, y0}, z);
dspxy = fnder(sp, [1, 1]);
fnplt(dspxy);
hold on;
syms x y;
z = (x^2-2*x)*exp(-x^2-y^2-x*y);
ezsurf(diff(diff(z, x), y), [-3, 3], [-2, 2])

样条插值的数值积分计算
数值积分运算
积分函数: $S_i$ = fnint(S)
积分函数的结果为样条对数
在区间[a, b]求取定积分
I = diff(fnval(fnint(S), [a, b]))

例8-18 正弦函数的数值积分
给定$\int_0^{\pi} \sin x dx$
生成一些稀疏的样本点
用样条积分的方程求出定积分及积分函数

x = [0, 0.4, 1 2, pi];
y = sin(x);
sp1 = csapi(x, y);
a = fnint(sp1, 1);
xx = fnval(a, [0, pi]);
I1 = xx(2) - xx(1);
sp2 = spapi(5, x, y);
b = fnint(sp2, 1);
xx = fnval(b, [0, pi]);
I2 = xx(2) - xx(1);
% 积分函数绘制
ezplot('-cos(t)+2', [0, pi]);
hold on;
fnplt(a, '--');
fnplt(b, ':');
hold off;

基于样条插值的数值积分运算
数值积分运算
积分函数 $S_i$ = fnint(S)
积分函数的结果为样条函数
在区间[a, b]求取定积分
I = diff(fnval(fnint(S), [a, b]))

例8-19 双重积分的积分曲面
积分公式 $J = \int_{-1}^{1} int_{-2}^{2} e^{-x^2/2} \sin(x^2+y) dx dy$
fnint()函数不能求解，可以考虑fnder()函数
微分阶次选为负值可以计算积分

x0 = -2: 0.1 : 2;
y0 = -1: 0.1: 1;
[x, y] = ndgrid(x0, y0);
z = exp(-x.^2/2).*sin(x.^2+y);
S = spapi({5, 5}, {x0, y0}, z);
S1 = fnder(S, [-1, -1]);
S2 = fnval(S1, {x0, y0});
surf(y0, x0, S2);

08.05 由已知数据拟合函数

多项式拟合
多项式拟合就是用多项式拟合已知数据
多项式拟合函数调用格式
p = polyfit(x, y, n)
其中：
x和y为原始的样本点构成的向量
n为选定的多项式阶次
p为多项式系数按降幂排列得出的行向量

例8-20 数据的多项式拟合
原型函数 $f(x) = (x^2 - 3x + 5)e^{-5x} \sin x$
使用该函数生成一些数据点
用多项式方法在不同的阶次下进行拟合
拟合该数据的3次多项式

x0 = 0: 0.1: 1;
y0 = (x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);
p3 = polyfit(x0, y0, 3);
x = 0: .01: 1;
ya = (x.^2-3*x+5).*exp(-5*x).*sin(x);
y1 = polyval(p3, x);
plot(x, y1, x, ya, x0, y0, 'o');

不同阶次拟合

p4 = polyfit(x0, y0, 4);
y2 = polyval(p4, x);
p5 = polyfit(x0, y0, 5);
y3 = polyval(p5, x);
p8 = polyfit(x0, y0, 8);
y4 = polyval(p8, x);
plot(x, ya, x0, y0, 'o', x, y2, x, y3, x, y4);

% Taylor幂级数展开
syms x;
y = (x^2-3*x+5)*exp(-5*x)*sin(x);
p = taylor(y, 'Order', 9)

例8-21 数据拟合
对下式进行多项式拟合，并观察拟合效果
$f(x) = 1/(1+25x^2), -1 \leq x \leq 1$
就不同的多项式阶次进行曲线拟合：

x0 = -1+2*[0:10]/10;
y0 = 1./(1+25*x0.^2);
x = -1: 0.01: 1;
ya = 1./(1+25*x.^2);
p3 = polyfit(x0, y0, 3);
y1 = polyval(p3, x);
p5 = polyfit(x0, y0, 5);
y2 = polyval(p5, x);
p8 = polyfit(x0, y0, 8);
y3 = polyval(p8, x);
p10 = polyfit(x0, y0, 10);
y4 = polyval(p10, x);
plot(x, ya, x, y1, x, y2, x, y3, x, y4);

函数线性组合的曲线拟合方法
已知函数的线性组合
$g(x) = c_1f_1(x) + c_2f_2(x) + c_3f_3(x) + \cdots + c_nf_n(x)$
Ac = y
$A = \begin{bmatrix}f_1(x_1) & f_2(x_1) & \cdots & f_m(x_1) \\ f_1(x_2) & f_2(x_2) & \cdots & f_m(x_2) \\ \vdots & \vdots & \ddots & \vdots \\ f_1(x_M) & f_2(x_M) & \vdots & f_m(x_M) \end{bmatrix}, y = \begin{bmatrix}y_1 \\ y_2 \\ \vdots \\ y_M \end{bmatrix}$

例8-21 多项式拟合
原型函数 $f(x) = (x^2 - 3x + t)e^{-5x} \sin x$
进行多项式拟合的各个函数为
$f_i(x) = x^{n+1-i}, i = 1, 2, \cdots, n$
观察多项式拟合效果

x = [0: .1: 2]';
n = 7;
A = [];
y = (x.^2+3*x+5).*exp(-5*x).*sin(x);
for i = 1:n+1
  A(:, i) = x.^(n+1-i);
end;
c = A\y

08.06 最小二乘曲线拟合

最小二乘曲线拟合
给定一组数据 $x_i, y_i, i = 1, 2, \cdots, N$
已知函数原型 $\hat y(x) = f(a, x), a 为待定系数$
最小二乘曲线拟合的目标为求出一组待定系数的值，使得目标函数为最小
$J = \min \limits _a \sum \limits_{i=1}^{N}[y_i - \hat y(x_i)]^2 = \min \limits_a \sum \limits_{i=1}^{N}[y_i - f(a, x_i)]^2$
matlab求解
[a, $j_m$] = lsqcurvefit(Fun, $a_0$, x, y, options)

例8-25 最小二乘曲线拟合
由下面的语句生成一组数据

x = 0: 0.1: 10;
y = 0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x);

该数据满足原型函数: $y(x) = a_1 e^{-a_2x} + a_3e^{-a_4x} \sin (a_5x)$

% 最小二乘拟合
f = @(a, x)a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x);
[xx, res] = lsqcurvefit(f, [1, 1, 1, 1, 1], x, y)

例8-26 最小二乘拟合
已知数据可能满足 $y(x) = ax + bx^2 e^{-cx} + d$
$\begin{array}{l|c c c c } \hline x_i & 0.1 & 0.2 & 0.3 &0.4 & 0.5 \\ y_i & 2.3201 & 2.6470 & 2.9707 & 3.2885 & 3.6008 \\ \hline x_i &0.6 & 0.7 & 0.8 & 0.9 & 1 \\ y_i & 3.990 & 4.2147 & 4.5191 & 4.8232 & 5.1275 \\ \hline \end{array}$
求满足数据的最小二乘解 a, b, c, d的值
令 $a_1 = a, a_2 = b, a_3 = c, a_4 = d$
原型函数 $y(x) = a_1x + a_2x^2e^{-a_3x} + a_4$

% 输入已知数据：
x = 0.1: 0.1: 1;
y = [2.3201, 2.6470, 2.9707, 3.2885, 3.6008, ...
  3.9090, 4.2147, 4.5191, 4.8232, 5.1275];
% 求解
f = @(a, x)a(1)*x + a(2)*x.^2 .* exp(-a(3)*x) + a(4);
a = [1;2;2;3];
[a, b, c, d] = lsqcurvefit(f, a, x, y);
a, b

多变量函数的最小二乘函数拟合
多变量待定函数 $z = f(z, X)$
目标函数：$z = f(a, x_1, x_2, \cdots, x_m)$
$J = \min \limits_{a} \sum \limits_{i=1}^{N} [y_i - \hat y(x_i)]^2 = \min \limits_{a} \sum \limits_{i=1}^{N} [y_i - f(a, x_i)]^2$
最小二乘拟合
$[a, j_m] = lsqcuryefit(Fun, a_0, X, y, options)$

例8-27 多变量函数拟合
原型函数 $z = a_1x^{a_2x} + a_3 y^{a_4(x+y)} + a_5z^{a_6(x+y+z)}$
数据文件:c8data1.dat
令 $x_1 = x, x_2 = y, x_3 = z$

f = @(a, X)a(1)*X(:, 1) .^ (a(2)*X(:, 1)) + ...
  a(3)*X(:, 2) .^ (a(4)*(X(:, 1) +X(:, 2))) + ...
  a(5)*X(:, 3) .^ (a(6)*(X(:, 1) + X(:, 2) + X(:, 3)));
% XX = load('E:\\Work\\201806\\math_4th\\c8data1.dat');
XX = load('c8data1.dat');
X = XX(:, 1:3);
v = XX(:, 4);
a0 = [2 3 2 1 2 3];
[a, fm] = lsqcurvefit(f, a0, X, v)

08.07 函数的有理式逼近

连分式的数学形式
$f(x) = b_1 + \dfrac{(x-a)^{c_1}}{b_2 + \dfrac{(x-a)^{c_2}}{b_3 + \dfrac{x-a)^{c_3}}{b_4 + \dfrac{(x-a)^{c_4}}{b_5 + \dfrac{(x-a)^{c_5}}{...}}}}}$
调用格式：cf = contfrac(f, n)
[cf, r] = contfrac(f, n, ‘x=a’)

function key=isanumber(a)
key=0; if length(a)~=1, return; end 
switch class(a)
    case 'double', key=1;
    case 'sym', try, double(a); key=1; catch, end
end

function [cf, r] = contfrac(f, varargin)
  [n, a] = default_vals({6,0}, varargin{:});
if isanumber(f)
    cf = feval(symengine, 'contfrac', f, n);
    p1 = char(cf);
    k = strfind(p1, ',');
    k1 = strfind(p1, '/');
    if nargout > 1
        r = sym(p1(k(end)+1:k1-1))/sym(p1(k1+1:end-1));
    end
else
    if isfinite(a)
        str=num2str(a);
    else
        str='infinity';
    end
    cf = feval(symengine, 'contfrac', f, ['x=' str], n);
    if nargout > 1
        r = feval(symengine, 'contfrac::rational', cf);
    end
end

例8-28 $\pi$的连分数近似
连分数近似cf = contfrac(pi, 20)
结果: 292大于其它值
截断于该项103993/33102
近似值: 3.141592653012

例8-29 函数的连分式近似
原型函数 $f(x) = \dfrac{e^{-x} \sin x}{(x+1)^3}$
10级连分式

syms x;
f = sin(x)*exp(-x)/(x+1)^3;
[cf, y] = contfrac(f, 10)
% 不同连分式近似
[cf1, f8] = contfrac(f, 8),
[cf2, f10] = contfrac(f, 10)

连分式拟合效果

% 拟合效果
ezplot(f, [0, 2]);
hold on;
ezplot(f8, [0, 2]);
ezplot(f10, [0, 2]);
figure;
ezplot(f, [0, 5]);
hold on;
ezplot(f8, [0, 5]);
ezplot(f10, [0, 5]);

% 参考点x = 0.5
[cf1, f8] = contfrac(f, 8, 0.5)
ezplot(f, [0, 5]);
hold on;
ezplot(f8, [0, 5]);
hold off;

有理式拟合：Pad$\acute{e}$近似
函数 f(s) 的幂级数展开
$f(s) = c_0 + c_1s + c_2s^2 + c_3s^3 + \cdots = \sum \limits_{i=0}^{\infty} c_i s^i$
$r/m阶的Pad\acute{e} 近似表示为：$
$G_m^{r}(x) = \dfrac{\beta_{r+1}s^r + \beta_{r} s^{r-1} + \cdots + \beta_1}{\alpha_{m+1} s^m + \alpha_m s^{m-1} + \cdots + \alpha_1} = \dfrac{\sum \limits_{i=1}^{r+1} \beta_i s^{i-1}}{\sum \limits_{i=1}^{m+1} \alpha_i s^{i-1}}$
$\alpha_1 = 1, \beta_1 = c_1$

系数计算
$系数 c_i 的递推计算 c_0 = b_{k+1},$
$c_i = b_{k+1-i} - \sum \limits_{j=0}^{i=1} c_j a_{n+1-i+j}, i = 1, 2, \cdots$
$设\sum \limits_{i=0}^{\infty} c_i s^i = G_m^r(s), 则可写出$
$\sum \limits_{i=1}^{m+1} \alpha_i s^{i-1} \sum \limits_{i=0}^{\infty} c_i s^{i} = \sum \limits_{i=1}^{r+1} \beta_i s^{i-1}$
令相应的 s 项系数的值相等，得出
$Wx = w, v = Vy$
其中:
$x = [\alpha_2, \alpha_3, \cdots, \alpha_{m+1}]^T$
$y = [\alpha_2, \alpha_3, \cdots, \alpha_{r+1}]^T$
$w = [-c_{r+2}, -c_{r+3}, \cdots, -c_{m+r+1}]^T$
$v = [\beta_2 - c_2, \beta_3 - c_3, \cdots, \beta_{r+1} - c_{r+1}]^T$
且有
$W = \begin{bmatrix}c_{r+1} & c_r & \cdots & 0 & \cdots & 0 \\ c_{r+2} & c_{r+1} & \cdots & c_1 & \cdots & 0 \\ \vdots & \vdots & \cdots & \vdots & \ddots & \vdots \\c_{r + m} & c_{r +m-1} & \cdots & c_{m-1} & \cdots & c_{r+1} \end{bmatrix}$

$V = \begin{bmatrix}c_1 & 0 & 0 & \cdots & 0 \\ c_2 & c_1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_r & v_{r-1} & c_{r-2} & \cdots & c_1 \end{bmatrix}$
编写函数计算
给定函数的$Pad\acute{e}$有理数近似调用格式
[n, d] = padefcn(c, r, m)

function [nP, dP] = padefcn(c, r, m)
  w = -c(r+2:m+r+1)';
  vv = [c(r+1:-1:1)';
  zeros(m-1-r,1)];
  W = rot90(hankel(c(m+r:-1:r+1),vv));
  V = rot90(hankel(c(r:-1:1)));
  x = [1 (W\w)'];
  y = [1 x(2:r+1)*V'+c(2:r+1)];
  dP = x(m+1:-1:1)/x(m+1);
  nP = y(r+1:-1:1)/x(m+1);

例8-30 函数的有理近似
给定函数 $f(x) = e^{-2x}$
求有理函数近似
分别选择分子阶次为0和不同的分母阶次

syms x;
c = taylor(exp(-2*x), 'Order', 10);
c = sym2poly(c);
c = c(end:-1:1);
ezplot(exp(-2*x), [0, 4]);
x = 0: 0.01: 8;
nd = [3:7];
for i = 1: length(nd)
  [n, d] = padefcn(c, 0, nd(i));
  y = polyval(n, x)./polyval(d, x);
  line(x, y);
end;

符号运算
符号运算版本的padefcnsym()函数

function c = polycoef(p, x)
  c=[];
  n = 0;
  p1=p;
  n1=1;
  nn=1;
  if nargin == 1
    x=symvar(p);
  end
  while (1)
    c = [c subs(p1, x, 0)];
    p1 = diff(p1, x);
    n = n+1;
    n1 = n1*n;
    nn = [nn, n1];
    if p1 == 0
      c=c./nn(1:end-1);
      c=c(end:-1:1);
      break;
    end
  end

function G = padefcnsym(f, r, m)
  c = taylor(f, 'Order', r+m+1);
  c = polycoef(c);
  c = c(end: -1: 1);
  w = -c(r+2:m+r+1)';
  vv = [c(r+1:-1:1)';
  zeros(m-1-r, 1)];
  W = rot90(hankel(c(m+r:-1:r+1), vv));
  V = rot90(hankel(c(r:-1:1)));
  X = [1 (W\w)'];
  y = [1 X(2:r+1)*V'+c(2:r+1)];
  dP = X(m+1:-1:1)/X(m+1);
  nP = y(r+1:-1:1)/X(m+1);
  syms x;
  G = poly2sym(nP, x)/poly2sym(dP, x);

例8-31 符号运算求取Pade近似
给定函数 $f(x) = e^{-2x}$
用符号运算方式求取 Pade近似模型

syms x;
f = exp(-2*x);
G = padefcnsym(f, 0, 8)

最终结果
$G(x) = \dfrac{315}{2(x^8+4x^7 + 4x^6 + 42x^5 + 105x^4 + 210x^3 + 315x^2 + 315x + 315/2)}$

08.08 特殊函数

误差函数与补误差函数
误差函数
$erf(z) = \dfrac{2}{\sqrt{\pi}}\int_0^{z} e^{-t^2} dt = \dfrac{1}{\sqrt{\pi}} \int_{-z}^{z} e^{-t^2} dt$
$erf(0) = \dfrac{2}{\sqrt{\pi}}\int_0^0 e^{-t^2} dt = 0, erf(\infty) = \dfrac{2}{\sqrt{\pi}} \int_0^{\infty} e^{-t^2} dt = 1, erf(-\infty) = -1$
补误差函数
$erfc(z) = \dfrac{2}{\sqrt{\pi}} \int_{z}^{\infty} e^{-t^2} dt$
两个函数的关系： erfc(z) = 1 - erf(z)

例8-34 误差函数与补误差函数的曲线
直接绘制：y = erf(z), y = erfc(z)

x = -5: 0.1: 5;
y1 = erf(x);
y2 = erfc(x);
plot(x, y1, x, y2, '--');

Gamma函数
Gamma函数的数学描述
$\Gamma(\alpha) = \int_0^{\infty} e^{-t}t^{\alpha - 1} d t$
Gamma函数调用格式: y = gamma(x)
其中，x可以为向量，矩阵或其它数据结构
gamma()函数可用于数据与符号运算
特殊性质：$\Gamma(n+1) = n!, \Gamma(\alpha + 1) = \alpha \Gamma(\alpha)$
常用性质：
$\Gamma(\dfrac{1}{2}) = \sqrt{\pi}, \Gamma(\alpha) \Gamma(1 - \alpha) = \dfrac{\alpha}{\sin \pi \alpha}$
$\Gamma(\alpha)\Gamma(-\alpha) = \dfrac{-\pi}{\alpha \sin \pi \alpha}$
$\Gamma(\dfrac{1}{2} + \alpha) \Gamma(\dfrac{1}{2} - \alpha) = \dfrac{\pi}{\cos \pi \alpha}$
$\lim \limits_{\alpha \to \infty} \dfrac{\alpha^z \Gamma(\alpha)}{\Gamma(\alpha+z)} = 1, Re z > 0$

例8-37 Gamma函数曲线
试绘制(-4, 3)区间内Gamma函数的曲线

x = -4: 0.0001: 3;
plot(x, gamma(x));
ylim([-40, 40]);

注意:
$\Gamma(\alpha+1) = \alpha \Gamma(\alpha), \Gamma(1) = 1$
$若\alpha为非负整数，则\Gamma(\alpha+1) = \alpha!$
$若\alpha为负整数, 则 \Gamma(\alpha+1) 趋于\pm \infty$

例8-38 不可积函数的解析求解
试写出下面无穷积分的“解析解”
$I_1 = \int_0^{\pi/2} \sin^{2m-1} t \cos ^{2n-1} t dt$
$I_2 = \int_0^{\infty} t^{x-1} cos t dt, x > 0$

syms m n t z;
syms x positive;
I1 = int(sin(t)^(2*m-1)*cos(t)^(2*n-1), t, 0, pi/2)
I2 = simplify(int(t^(x-1)*cos(t), t, 0, inf))

不完整Gamma函数
不完整$\Gamma$-函数的数学描述
$\Gamma(x, \alpha) = \dfrac{1}{\Gamma(\alpha)} \int_0^{x} e^{-t} t^{\alpha-1} dt, x \geq 0$
求取不完整$\Gamma$-函数的调用格式
y = gammainc(x, b)

Beta函数
Beta函数的数学描述
$B(z, m) = \int_0^1 t^{z-1} (1 - t)^{m-1} dt, Re(m) > 0, Re(z) > 0$
$B(z, m) = \dfrac{\Gamma(z) \Gamma(m)}{\Gamma(z + m)}$
Beta函数的调用格式： y = beta(m, z)

Bessel函数
对于Bessel微分方程
$t^2 \dfrac{d^2 x}{dt^2} + t \dfrac{dx}{dt} + (t^2 - \lambda^2) x = 0$
该方程的通解为
$x(t) = C_1J_{\lambda}(t) + C_2J_{-\lambda}(t)$
$其中, C_1和C_2为任意常数$

Bessel函数计算
第一类Bessel函数调用格式: y = besselj($\lambda, x$)
第二类Bessel函数调用格式: y = bessely($\lambda, x$)
第三类Bessel函数调用格式: y = besselh($\lambda, x$)
其中, $\lambda$为阶次

Legendre函数
对于Legendre微分方程
$(1-t^2)\dfrac{d^2 x}{dt^2} - 2t \dfrac{dx}{dt} + n(n+1) x = 0$
该方程的通解为
$x(t) = C_1P_n(t) + C_2Q_n(t)$
$其中, C_1和C_2为任意常数$

Legendre函数计算
计算关联Legendre函数 $P_n^m(t)$ 调用格式: Y = legendre(n, x)
其中, Y为一个矩阵，其各行分别为
$P_n^0(x), P_n^1(x), \cdots, P_n^n(x)$
该函数要求：$-1 < x < 1$

例8-41 Legendre 函数曲线

x = -1: 0.04: 1;
Y = legendre(2, x);
plot(x, Y);

超几何函数
超几何函数的一般形式
$_pF_q(a_1, \cdots, a_p; b_1, \cdots, b_q; z) = \dfrac{\Gamma(b_1)\cdots \Gamma(b_q)}{\Gamma(a_1)\cdots \Gamma(a_p)} \sum \limits_{k=0}^{\infty} \dfrac{\Gamma(a_1+k)\cdots \Gamma(a_p+k) z^k}{\Gamma(b_1+k)\cdots \Gamma(b_q+k) k!}$
符号运算: y = hypergeom($[a_1, \cdots, a_p], [b_1, \cdots, b_q], z$)
数值运算: $y_1$ = subs($y, z, x_0$)

08.09 Mittag-Leffler函数

Mittag-Leffler函数
指数函数的幂级数展开
$e^z = \sum \limits_{k=0}^{\infty} \dfrac{z^k}{k!} = \sum \limits_{k=0}^{\infty} \dfrac{z^k}{\Gamma(k+1)}$
单参数 Mittag-Leffler 函数 $E_1(z) = e^z$
$E_{\alpha}(z) = \sum \limits_{k=0}^{\infty} \dfrac{z^k}{\Gamma(\alpha k + 1)}, Re(\alpha) > 0$
其它常用的单参数 Mittag-Leffler 函数
$E_2(z) = \sum \limits_{k=0}^{\infty} \dfrac{z^k}{\Gamma(2k+1)} = \sum \limits_{k=0}^{\infty} \dfrac{(\sqrt{z})^{2k}}{(2k)!} = \cosh \sqrt{z}$
$E_{1/2}(z) = \sum \limits_{k=0}^{\infty} \dfrac{z^k}{\Gamma(k/2+1)} = e^{z^2}(1 + erf(z)) = e^{z^2} erfc(-z)$

双参数 Mittag-Leffler 函数 $Re(\alpha) > 0, Re(\beta) > 0$
$E_{\alpha, \beta}(z) = \sum \limits_{k=0}^{\infty} \dfrac{z^k}{\Gamma(\alpha k +\beta)}, E_{\alpha, 1}(z) = E_{\alpha}(z)$

双参数Mittag-Leffler函数的一些公式
$E_{1, 2}(z) = \sum \limits_{k=0}^{\infty} \dfrac{z^k}{\Gamma(k+2)} = \dfrac{1}{z} \sum \limits_{k=0}^{\infty} \dfrac{z^{k+1}}{(k+1)!} = \dfrac{e^z - 1}{z}$
$E_{1, 3}(z) = \sum \limits_{k=0}^{\infty}\dfrac{z^k}{\Gamma(k+3)} = \dfrac{e^z - 1 - z}{z^2}$
$E_{1, m}(z) = \sum \limits_{k=0}^{\infty} \dfrac{z^k}{\Gamma(k+m)} = \dfrac{1}{z^{m-1}}(e^z - \sum \limits_{k=0}^{m-2}\dfrac{z^k}{k!})$
$E_{2,2}(z) = \sum \limits_{k=0}^{\infty}\dfrac{z^k}{\Gamma(2k+2)} = \dfrac{1}{z} \sum \limits_{k=0}^{\infty} \dfrac{(\sqrt{z})^{2k+1}}{(2k+1)!} = \dfrac{\sinh \sqrt{z}}{\sqrt{z}}$

% Mittag-Leffler解析函数
function f = mittag_leffler(aa, z)
  aa = [aa 1];
  a = aa(1);
  b = aa(2);
  syms k;
  f = simplify(symsum(z^k/gamma(a*k+b), k, 0, inf));

函数入口:
f1 = mittag_leffler($\alpha$, z)
f1 = mittag_leffler([$\alpha, \beta$], z)

例8-36 Mittag-Leffler函数
单参数 Mittag-Leffler 函数推导
$\alpha = 1/3, 3, 4, 5$

syms z;
I1 = mittag_leffler(1/sym(3), z),
I2 = mittag_leffler(3, z)
I3 = mittag_leffler(4, z)
I4 = mittag_leffler(5, z)
% 在matlab2008a下可以求出非超几何函数

例8-37 双参数 Mittag-Leffler 函数推导

syms z;
I5 = mittag_leffler([4, 1], z),
I6 = mittag_leffler([4, 5], z),
I7 = mittag_leffler([5, 6], z),
I8 = mittag_leffler([1/sym(2), 4], z)

双参数 Mittag-Leffler 函数的整数阶导数
$\dfrac{d^n}{d z^n} E_{\alpha, \beta}(z) = \sum \limits_{k=0}^{\infty} \dfrac{(k+n)!}{k!\Gamma(\alpha k + \alpha n + \beta)} z^k$

调用格式:
f = ml_func($\alpha, z, n, \epsilon_0$)
f = ml_func($[\alpha, \beta], z, n, \epsilon_0$)

function f = ml_func(aa, z, varargin)
aa=[aa,1,1,1]; a=aa(1); b=aa(2); c=aa(3); q=aa(4); f=0; k=0; fa=1; 
[n,eps0]=default_vals({0,eps},varargin{:});
if n==0
   while norm(fa,1)>=eps0
      fa=gamma(k*q+c)/gamma(c)/gamma(k+1)/gamma(a*k+b) *z.^k;
      f=f+fa; k=k+1;
   end
   if ~isfinite(f(1))
      if c*q==1
         f=mlf(a,b,z,round(-log10(eps0))); f=reshape(f,size(z));
      else, error('Error: truncation method failed'); end, end
else
   aa(2)=aa(2)+n*aa(1); aa(3)=aa(3)+aa(4)*n;
   f=gamma(q*n+c)/gamma(c)*ml_func(aa,z,0,eps0);
end

例8-38 Mittag-Leffler函数
得出 $E_{1/2, 4}(z)$ 的解析解与数值解(比较耗时)

syms z;
I8 = mittag_leffler([1/sym(2), 4], z);
t = 0: 0.01: 2;
y = subs(I8, z, t);
y1 = ml_func([1/2, 4], t);
plot(t, y, t, y1, '--');

% 数值解
t = 0: 0.01: 0.5;
y1 = ml_func([1/2, 4], t);
plot(t, y1);

例8-46 绘制函数 $E_{0.8, 0.9}(z)的实部与虚部表面图$
假设 z 是复数变量

% 生成网格, 构造复数矩阵, 绘制实部表面图
[x, y] = meshgrid(-6: 0.2: 6);
z = x + sqrt(-1)*y;
L = ml_func([0.8, 0.9], z);
L1 = real(L);
ii = find(L1 > 3);
L1(ii) = 3;
ii = find(L1 < -3);
L1(ii) = -3;
surf(x, y, L1);
axis([-6 6 -6 6 -3 3])
% 虚部表面图
figure;
L2 = imag(L);
ii = find(L2 > 3);
L2(ii) = 3;
ii = find(L2 < -3);
L2(ii) = -3;
surf(x, y, L2);
axis([-6 6 -6 6 -3 3]);

多参数Mittag-Leffler函数
三参数Mittag_leffler函数
$E_{\alpha, \beta}^{\gamma}(z) = \sum \limits_{k=0}^{\infty} \dfrac{(\gamma)_k}{\Gamma(\alpha k + \beta)} \dfrac{z^k}{k!}, Re(\cdot) > 0$
四参数Mittag_leffler函数
$E_{\alpha, \beta}^{\gamma, q}(z) = \sum \limits_{k=0}^{\infty} \dfrac{(\gamma)^{kq}}{\Gamma(\alpha k + \beta)} \dfrac{z^k}{k!}, q \in N$
Pochhammer符号
$(\gamma)_k = \gamma(\gamma+1)(\gamma+2)\cdots(\gamma+k-1) = \dfrac{\Gamma(k + \gamma)}{\Gamma(\gamma)}$
$E_{\alpha, \beta}(z) = E_{\alpha, \beta}^{1}(z), E_{\alpha, \beta}^{\gamma}(z) = E_{\alpha, \beta}^{\gamma, 1}(z)$

Mittag-Leffler函数的导数
n阶导数公式
$\dfrac{d^n}{dz^n}E_{\alpha, \beta}^{\gamma, q}(z) = (\gamma)_{qn} E_{\alpha, \beta+n\alpha}^{\gamma + qn, q}(z)$
matlab求解：
f = ml_func($\alpha, z, n, \epsilon_0$)
f = ml_func($[\alpha, \beta], z, n, \epsilon_0$)
f = ml_func($[\alpha, \beta, \gamma], z, n, \epsilon_0$)
f = ml_func($[\alpha, \beta, \gamma, q], z, n, \epsilon_0$)

08.10 信号的相关分析

信号的相关分析
信号$x(t)$的自相关函数的定义：
$R_{xx}(\tau) = \lim \limits_{T \to \infty} \dfrac{1}{T} \int_0^{T} x(t) x(t + \tau) d t, \tau \geq 0$
且相关函数为偶函数，即 $R_{xx}(-\tau) = R_{xx}(\tau)$
两个信号$x(\tau) 和 y(\tau)$ 的互相关函数为：
$R_{xy}(\tau) = \lim \limits_{T \to \infty} \dfrac{1}{T} \int_0^{T} x(t) y(t + \tau) dt, \tau \geq 0$

例8-39 自相关函数计算
给定信号$x(t) = A_1 \cos(\omega_1 t + \theta_1) + A_2 \cos(\omega_2 t + \theta_2)$
求它的自相关函数

syms A1 A2 w1 w2 t1 t2 t tau T;
x = A1*cos(w1*t + t1) + A2*cos(w2*t + t2);
Rxx = simplify(limit(int(x*subs(x, t, t+tau), t, -T, T)/2/T, T, inf))

结果:
$R_{xx}(\tau) = \dfrac{1}{2}A_1^2 \cos(\omega_1 \tau) + \dfrac{1}{2} A_2^2 \cos(\omega_2 \tau)$

相关分析的数值方法
在实验中测出两组数值: $x_i, y_i, (i = 1, 2, \cdots, n)$
计算两组数据的相关系数
$\gamma_{xy} = \dfrac{\sqrt{\sum (x_i - \bar x)(y_i - \bar y)}}{\sqrt{(x_i - \bar x)^2} \sqrt{(y_i - \bar y)^2}}$
求取已知向量的相关系数矩阵
R = corrcoef(x, y)
R = corrcoef([x, y])

例 8-40 相关系数矩阵
用下列两个原型函数分别生成一组数据, 并求取其相关系数矩阵
$y_2 = t e^{-4t} \cos 3t 和 y_1 = t e^{-4t} \sin 3t$

x = 0: 0.01: 5;
y1 = x.* exp(-4*x) .* sin(3*x);
y2 = x .* exp(-4*x) .* cos(3*x);
R = corrcoef(y1, y2)

相关函数计算
在实验中测出两组数据 $x_i, y_i, (i = 1, 2, \cdots, n)$
定义 x 序列的自相关函数
$c_{xx}(k) = \dfrac{1}{N} \sum \limits_{l=1}^{n-[k]-1} x(l)x(k+l), 0 \leq k \leq m-1$
其中, m < n, 自相关函数是偶函数
定义出相关函数
$c_{xy}(k) = \dfrac{1}{N}\sum \limits_{l=1}^{n-[k]-1} x(l) y(k+l), 0 \leq k \leq m-1$

求自相关的调用格式
$C_{xx}$ = xcorr(x, N)
求互相关的调用格式
$C_{xx}$ = xcorr(x, y, N)
其中, N为k的最大取值, 可以忽略

例8.41 相关函数
用函数$y_1 = te^{-4t}\sin 3t, y_2 = te^{-4t} \cos 3t$
生成一组数据，用数值方法计算相关函数
并和已知理论曲线进行比较

t = 0: 0.01: 5;
x = t.*exp(-4*t).*sin(3*t);
y = t.*exp(-4*t).*cos(3*t);
N = 150;
c1 = xcorr(x, N);
x1 = [-N: N];
stem(x1, c1);
figure;
c1 = xcorr(x, y, N);
x1 = [-N: N];
stem(x1, c1)

08.11 信号滤波与滤波器设计

例8-42 污染的信号
给定信号 $y(x) = e^{-x} \sin (5x)$
叠加标准误差为 $\sigma = 0.05$ 的零均值的白噪声信号
绘制噪声污染后的信号曲线

x = 0: 0.002: 2;
y = exp(-x).*sin(5*x);
r = 0.05*randn(size(x));
y1 = y + r;
plot(x, y1);

线性滤波器的一般模型
离散线性时不变滤波器模型一般表达式为:
$H(z) = \dfrac{b_1 + b_2 z^{1} + b_3 z^{-2} + \cdots + b_{n+1} z^{-n}}{1 + a_1z^{-1} + a_2z^{-2} + \cdots + a_m z^{-m}}$
假设输入信号为 x(n), 则经过滤波器后的输出信号为：
$y(k) = -a_1 y(k-1) - \cdots - a_my(k-m) + b_1x(k) + b_2x(k-1) + \cdots + b_{n+1}x(k-n)$

三种滤波器类型
根据n和m的不同取值，定义常用的滤波器
FIR滤波器：又称为有限长脉冲响应滤波器, m = 0, a为标量，称为移动平均模型(MA)
IIR滤波器：又称全极点无限长脉冲响应滤波器，也称自回归模型(AR)，n = 0, b为标量
ARMA滤波器：又称为一般IIR滤波器，也称为自回归移动平均模型(ARMA),要求n与m均不为0

滤波的信号向量调用格式: y = filter(b, a, x)
a,b 向量表示滤波器， x为需要过滤的信号
对滤波器进行放大倍数分析
[h, w] = freqz(b, a, N)
只获得放大倍数的幅值
plot(w, abs(h))

例8-43 滤波器效果
假设已经设计滤波器H(z)如下, 试得出该滤波器的放大倍数, 并观察滤波后的信号
$A(z) = 1.2296 \times 10^{-6}(1 + z^{-1})^7$
$B(z) = 1- 0.7265z^{-1}$
$C(z) = 1 - 1.488z^{-1} + 0.5641z^{-2}$
$D(z) = 1 - 1.595z^{-1} + 0.6769z^{-2}$
$E(z) = 1 - 1.78 z^{-1} + 0.8713z^{-2}$
$H(z) = \dfrac{A(z)}{B(z)C(z)D(z)E(z)}$

b = 1.2296e-6*conv([1 4 6 4 1], [1 3 3 1]);
a = conv([1, -0.7265], ...
   conv([1, -1.488, 0.5644], ...
   conv([1, -1.595, 0.6769], ...
   [1, -1.78, 0.8713])));
x = 0: 0.002: 2;
y = exp(-x).*sin(5*x);
r = 0.05*randn(size(x));
y1 = y+r;
[h, w] = freqz(b, a, 100);
plot(w, abs(h))
figure;
y2 = filter(b, a, y1);
plot(x, y1, x, y2);

滤波器及matlab实现
常用的滤波器类型
Butterworth滤波器
[b, a] = butter(n, $w_n$)
Chebyshev I 型滤波器
[b, a] = cheby1(n, r, $w_n$)
Chebyshev II型滤波器
[b, a] = cheby2(n, r, $w_n$)

例8-44 滤波器仿真
给定信号 $y(x) = e^{-x} \sin (5x)$
叠加标准差为 $\sigma = 0.05$ 的零均值的白噪声信号
试对阶次及不同的$w_n$ 值, 设计Butterworth 滤波器并比较滤波效果 $w_n = 1$

f1 = figure;
f2 = figure;
for n = 5: 2: 20
  figure(f1);
  [b, a] = butter(n, 0.1);
  y2 = filter(b, a, y1);
  plot(x, y2);
  hold on;
  figure(f2);
  [h, w] = freqz(b, a, 100);
  plot(w, abs(h));
  hold on;
end;

阶次为7时

f1 = figure;
f2 = figure;
for wn = 0.1: 0.1: 0.7
    figure(f1);
    [b, a] = butter(7, wn);
    y2 = filter(b, a, y1);
    plot(x, y2);
    hold on;
    figure(f2);
    [h, w] = freqz(b, a, 100);
    plot(w, abs(h));
    hold on;
end

其它滤波器
设计高通滤波器
[b, a] = butter(n, $w_n$, ‘high’)
设计带通滤波器
[b, a] = butter(n, [$w_1, w_2$])