函数式编程——巧用匿名函数在m程序中返回复杂表达式的函数句柄

1 引言

matlab的优点在于利用矩阵完成数值计算，但在数学分析中，我们需要的不仅仅是一些数，还需要一个解析的函数表达式，来单独对某些变量进行映射分析。而且，matlab中内置了大量的科研工具函数，如何充分利用前人写好的函数来快速生成自己想要的函数，以实现更好的调用？匿名函数的出现为上述问题提供了很好的解决方法。

2 问题背景

下面以矩形波函数展开成傅里叶级数的例题，来说明几种比较不错的使用方式。
在高数第六版下册（同济大学）的321页例3中，把宽为τ、高为h，周期为T的矩形波（图12-18）展开成了复数形式的傅里叶级数，即把（下面只写出一个周期的表达式）
$$u(t)=\{\begin{array}{c}0,-\frac{T}{2}\le t<-\frac{\tau }{2}\\ h,-\frac{\tau }{2}\le t<\frac{\tau }{2}\\ 0,\frac{\tau }{2}\le t<\frac{T}{2}\end{array}$$
变为
$$U\left(t\right)=\frac{h\tau }{T}+\frac{h}{\pi }\sum _{\begin{array}{c}n=-\infty \\ n\ne 0\end{array}}^{\infty }\frac{1}{n}\sin \frac{n\pi \tau }{T}{e}^{i\frac{2n\pi t}{T}}$$
那么，如何在matlab中得到这两个表达形式不同的解析函数？

3 原始函数

3.1 利用已有函数

在matlab的命令窗口中输入edit square，回车可以看到内置方波函数square(t,duty)的源代码。

function s = square(t,duty)
if nargin < 2
    duty = 50;
end
% 't' input must be real
validateattributes(t,{'numeric'},{'real'},mfilename,'t',1);
% 'duty' input must be a real scalar
validateattributes(duty,{'numeric'},{'real','scalar'},mfilename,'duty',2);
% Compute values of t normalized to (0,2*pi)
tmp = mod(t,2*pi);
% Compute normalized frequency for breaking up the interval (0,2*pi)
w0 = 2*pi*duty(1,1)/100;
% Assign 1 values to normalized t between (0,w0), 0 elsewhere
nodd = (tmp < w0);
% The actual square wave computation
s = 2*nodd-1;

可以用如下调用实例来作图感受一下（为方便起见，以下作图都只展示矩形波的3个波峰）。

t = -.05:.0001:.0525;
y = square(2*pi*30*t);plot(t,y);
set(gca,'YLim',[-1.2 1.2]);

进一步分析可以知道，函数square是将如下函数f(t)在整个实轴上进行了周期延拓。
$$f(t)=\{\begin{array}{c}1,0\le t<2\pi \times duty\\ -1,2\pi \times duty\le t<2\pi \end{array}$$

既然有了源码，我们完全可以将square函数修改后另存为符合例题要求的原函数，但本文将采取另一种方式，即只对square函数进行一下“外部包装”。特别是当某些Built-in函数不能查看源码时，这种方式更能派上用场。下面先上写好的orig.m文件：

function u= orig(h,k,T)
%h,k,T分别等于例题中的h,τ,T
u=@(t)square(2*pi/T*(t+k/2),k/T*100);
u=@(t)(u(t)+1)/2*h;
end

需要注意的是，orig函数的输出变量u并非数值，而是一个函数句柄！也就是说，该m文件的作用是：在边界条件h,k,T给定的情况下，用orig函数将@square映射为函数句柄u！
orig函数中第一句是将square的函数图像左移τ/2，以便和书中的矩形波一样关于竖轴对称，然后再横向进行系数为2π/T的伸缩变形，同时在square函数中设置占空比为k/T*100%；orig函数中第二句是将square的函数图像上移1个单位，以便和书中的矩形波一样落在上半平面，然后再竖向进行系数为h/2的伸缩变形。
于是，有了orig.m文件后，我们便可以利用它输出的符合例题要求的u(t)了。例如，宽为0.4、高为1，周期为1的矩形波可以这样作出：

u=orig(1,0.4,1);
t=-1.5:0.001:1.5;plot(t,u(t));set(gca,'YLim',[0 1.2])

3.2 分段函数

其实读完square函数的源码后，不难发现其核心算法——将时间t除以周期的余数与占空时间进行比较。借鉴这一点，我们完全可以直接将函数u(t)写成分段形式。下面先上写好的origd.m文件：

function u = origd(h,k,T)
%h,k,T分别等于例题中的h,τ,T
u=@(t)(mod(t,T)<k)*1+(mod(t,T)>=k)*0;
u=@(t)u(t+k/2);
end

需要注意的是，origd函数中第一句利用了比较大小的返回值——若为真返回1，若为假返回0。而且，如果比较语句括号外紧跟的是一个含t表达式，则 * 应该改成 .*，不然当t的输入值为矩阵或向量时，会产生一个典型的错误：

Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns
in the first matrix matches the number of rows in the second matrix. To perform
elementwise multiplication, use '.*'.

显然，这个函数句柄u与上一节中的在数学上是等价的。

4 傅里叶级数

4.1 函数的累计

对于含求和符号的傅里叶级数表达式，如果仅仅是计算出数值，在不考虑快速傅里叶变换时，毫无疑问采用行向量乘以列向量的形式更有效率，即将表达式看成
$$U\left(t\right)=\frac{h\tau }{T}+\frac{h}{\pi }\left[\begin{array}{ccc}\begin{array}{cc}\cdots & -\sin \frac{-\pi \tau }{T}\end{array}& \sin \frac{\pi \tau }{T}& \cdots \end{array}\right]\left[\begin{array}{c}\begin{array}{c}⋮\\ e^{i\frac{-2\pi t}{T}}\end{array}\\ \begin{array}{c}{e}^{i\frac{2\pi t}{T}}\\ ⋮\end{array}\end{array}\right]$$

但为了得到解析函数句柄，我们不得不牺牲matlab的矩阵优势，采用for循环对匿名函数进行累加操作。

function [u,c0,c] = rect(h,k,T,n,t)
%h,k,T,n,t分别等于例题中的h,τ,T,n,t
%当输入参数小于5个,即不输入t时，u为傅里叶级数的函数句柄，否则u为利用傅里叶级数算出的时域数值
%c0为直流分量,c为复数形式傅里叶级数中各指数项系数
a=[-n:-1,1:n];
c0=h*k/T;
c=h/pi./a.*sin(a*pi*k/T);
if nargin<5
    u=@(t)c0;         %用恒定值c0初始化函数u(t)
    for j=1:n
        u=@(t)u(t)+c(n+j)*exp(1i*2*j*pi*t/T)+c(n+1-j)*exp(-1i*2*j*pi*t/T);
    end
else
    u=c0+c*exp(1i*2*pi/T*a'*t);
end    
end

上述rect.m中整合了两种计算方式：当不输入时间向量数值t时，返回的u为傅里叶级数的函数句柄；当输入时间向量数值t时，返回的u为对应的时域信号值。即在数学计算上，前者先采用for循环对函数U(t)进行累加再代入t值进行求解，后者先直接代入t值再采用矩阵乘法计算出U(t)。也就是说，在命令窗口中输入

[u,c0,c]=rect(1,0.4,1,800);t=-1.5:0.0001:1.5;plot(t,u(t));

和

t=-1.5:0.0001:1.5;[u,c0,c]=rect(1,0.4,1,800,t);plot(t,u);

是等价的，两者均可以得到下图。
有趣的是，后者会比前者多输出一个warning:

Imaginary parts of complex X and/or Y arguments ignored.

造成这个warning的原因，很可能是matlab在对复杂函数U(t)调用前对表达式进行了整合简化，消去了复数形式的傅里叶级数中虚数i，从而在实际计算时，采用的是实数形式的傅里叶级数表达式
$$U\left(t\right)=\frac{a_0}{2}+\sum _{n=1}^{\infty }\left(a_n\cos \frac{n\pi t}{L}+b_n\sin \frac{n\pi t}{L}\right)$$
其中L=T/2。
而利用rect函数直接计算出对应时间的U(t)时，每一步都严格按照复数形式进行数值计算，从而使每一个结果都多出了很小的虚部，通过查看imag(u)的值可以知道，该虚部均保持在10^-16数量级。
但即便如此，在matlab中后者的运行效率依然领先于前者，分别在两者的前面添加语句“tic”、后面添加语句“toc”，可以得到Elapsed time分别为8.9s和0.78s。如果再为后者加上快速傅里叶变换的翅膀，两者就更不可同日而语了。这反过来正验证了一句经验之谈：
要想提高matlab中程序的执行效率，那么越是大型的数值计算，越是要用矩阵运算代替循环！