运用Matlab验证卷积并引出傅里叶变换

冲浪的时候看到一个题目，发现并不符合先前做题的规律：

如图，按以往的做题经验，两个函数左，右坐标之和应该等于卷积之后的左右坐标。但是这个函数并不符合。

后来发现其实只是因为这个未知的函数卷积之后使得（0，2）以外的区间的值全部抵消接近于零。

接下来将用matlab演示和验证

首先先定义自变量区间,matlab具体语法这里不再赘述

interval相当于微分的间距，取0.001为间距进行微分。

至于其他相关参数我直接取距离为50，20（数据太大计算机算力不够就会算的很慢）

T0=50;
intervalT = 0.001;
t = 0 : intervalT : T0;

W0 = 20;
intervalW=0.001;
w = 0: intervalW : W0;

创建函数function(W)

FunctionW = zeros(1,length(w));

这里先分析前半个函数，讨论该函数无意义的点。用一个循环来筛选

for i=1:length(w)
    if sin(w(i))>1e-6%判断分母是否为0导致函数不存在
        functionW(i)=sin(0.5*w(i))^2/(w(i).*sin(w(i)));%".*"为矩阵乘法
    else
        functionW(i)=0.25;%趋于零的点取极限值为0.25
    end
end

注意function(w)不带cos，单纯前半部分。

很容易得到，上图中两个函数都是偶函数，所以h(t)也是偶函数。所以我们可以改变积分的上下限以达到近似处理的目的。（这种处理会对图像造成误差，但是无法避免，因为电脑无法处理无限）。

接下来再用嵌套循环求h(t)

%求h(t)函数
h =zeros(1,length(t));
for i=1:length(t)
    for j=1:length(functionW)
        h(i)=h(i)+functionW(j)*cos(w(j)*t(i))*intervalW;
    end
end
h=2*h/pi%把上面的结果乘以pai分之2

求出正半轴的值，运用偶函数的性质即可得到负半轴

%由于h(t)为偶函数，所以只需要将函数翻转到负半轴即可
t=[-t(end:-1:2),t];%t的取值由-t到t，end指t最后的取值，从最后的取值，每次跨度-1，一直到2。
h=[h(end:-1:2),h];%h(t)为偶函数，所以无负号，其余一样。

画图

subplot(2,1,1);
%subplot是将多个图画到一个平面上的工具。其中，m表示是图排成m行，n表示图排成n列，
%也就是整个figure中有n个图是排成一行的，一共m行，如果m=2就是表示2行图。
%p表示图所在的位置，p=1表示从左到右从上到下的第一个位置。

plot(t,h);%直接将两个坐标画出来

验证是否能正推逆推，利用卷积公式求y(t)

%验证h(t)和x(t)卷积之后得到的y(t)是什么样子
newT=-10:intervalT:10;
y=zeros(1,length(newT));
for i=1:length(newT)
    time=newT(i);
    for tao=-20:0.001:20
        hindex=1+round((tao-t(1))/intervalT);
        if time-tao > 0 && time-tao < 2 && hindex >= 1 && hindex <= length(h)
            y(i) = y(i)+h(hindex)*0.001;
        end
    end
end

最后画出图像，证毕

%画y(t)
subplot(2,1,2);
plot(newT,y);
    

图像生成结果：

h(t)

y(t) 因为用了三次近似，所以会有误差，但是大致形状一致。

可以看到在0~2波形明显，其余区间几乎为0。

完整代码：

T0=50;
intervalT = 0.001;
t = 0 : intervalT : T0;

W0 = 20;
intervalW=0.001;
w = 0: intervalW : W0;

FunctionW = zeros(1,length(w));

for i=1:length(w)
    if sin(w(i))>1e-6%判断分母是否为0导致函数不存在
        functionW(i)=sin(0.5*w(i))^2/(w(i).*sin(w(i)));%".*"为矩阵乘法
    else
        functionW(i)=0.25;%趋于零的点取极限值为0.25
    end
end

%求h(t)函数
h =zeros(1,length(t));
for i=1:length(t)
    for j=1:length(functionW)
        h(i)=h(i)+functionW(j)*cos(w(j)*t(i))*intervalW;
    end
end
h=2*h/pi;
%由于h(t)为偶函数，所以只需要将函数翻转到负半轴即可
t=[-t(end:-1:2),t];%t的取值由-t到t，end指t最后的取值，从最后的取值，每次跨度-1，一直到2。
h=[h(end:-1:2),h];%h(t)为偶函数，所以无负号，其余一样。

%画图
subplot(2,1,1);
%subplot是将多个图画到一个平面上的工具。其中，m表示是图排成m行，n表示图排成n列，
%也就是整个figure中有n个图是排成一行的，一共m行，如果m=2就是表示2行图。
%p表示图所在的位置，p=1表示从左到右从上到下的第一个位置。

plot(t,h);%直接将两个坐标画出来

%验证h(t)和x(t)卷积之后得到的y(t)是什么样子
newT=-10:intervalT:10;
y=zeros(1,length(newT));
for i=1:length(newT)
    time=newT(i);
    for tao=-20:0.001:20
        hindex=1+round((tao-t(1))/intervalT);
        if time-tao > 0 && time-tao < 2 && hindex >= 1 && hindex <= length(h)
            y(i) = y(i)+h(hindex)*0.001;
        end
    end
end
%画y(t)
subplot(2,1,2);
plot(newT,y);