傅里叶级数3D动画演示，编写环境：MATLAB2021b
classdef Msg < event.EventData
%MSG 定义事件消息，event.EventData子类，消息中封装有要发布的数据
%Data 数据可以是任意类型
properties
Data %没有定义类型，所以发布的数据可以是任意类型
end
methods
function obj = Msg(newData)
%MSG 事件消息构造函数，由传递的数据决定数据类型
obj.Data = newData;
end
end
end
classdef Prt < handle
%PRT 信号参数类
% 作者：易
% 964615453@qq.com
% 2022-5-7
properties
SingleName %信号名称
T %信号周期
N %谐波次数
t %采样时间序列
y %信号函数表达式句柄，例如：y = @(t)square(w.*t)
SR='Seconds' %时间坐标单位
enRadio = false %是否出发广播帧，ture存储GIF动画文件
nD='view(app.axh,3)' %'view(app.axh,3)'绘制3维图动画，'view(app.axh,2)'绘制2维动画
end
methods
function P = Prt(SingleName,T,N,t,y,SR,enRadio,nD)
%构造函数，并参数赋值
P.SingleName=SingleName;
P.T=T;
P.t=t;
P.N=N;
P.y=y;
P.SR=SR;
P.enRadio=enRadio;
P.nD=nD;
end
end
end
classdef GenerateSingls < Prt
%GENERATESINGL 常见信号生成
%参数含义见Prt.m
% 作者：易
% 964615453@qq.com
% 使用枚举不是最好的选择，之前没有用过枚举所以尝试使用体会一下
% 2022-5-11
% 2022-5-15 为了演示旋转动画效果增加了3个旋转的方波
enumeration
Sawtooth('锯齿波',4,10,-4:0.025:4,@(t)sawtooth(2*pi/4.*t,.745),'Seconds',false,'view(app.axh,[30,35])')
Square('方波',4,20,-4:0.025:4,@(t)square(2*pi/4.*t),'Seconds',false,'view(app.axh,[30,35])')
%Square30('方波旋转30',4,20,-4:0.025:4,@(t)(exp(1j*pi/6).*square(2*pi/4.*t)),'Seconds',false,'view(app.axh,[30,35])')
%Square45('方波旋转45',4,20,-4:0.025:4,@(t)(exp(1j*pi/4).*square(2*pi/4.*t)),'Seconds',false,'view(app.axh,[30,35])')
%Square90('方波旋转90',4,20,-4:0.025:4,@(t)(exp(1j*pi/2).*square(2*pi/4.*t)),'Seconds',false,'view(app.axh,[30,35])')
%================================================================
%作者注：
%以上3条方波的语句可以取消注释，重启Main，选择其中的信号并启动3D动画
% 转动3D图到到旋转的角度，试试看；
% 有点烧脑，由于投影的切面始终是正交方向，所以看着有点奇怪。
%这应该如何理解呢？电磁波的极化角吗？
%===============================================================
Triangle('三角波',4,10,-4:0.025:4,@(t)mod(t,4),'Seconds',false,'view(app.axh,[30,35])')
AbsSin('正正弦波',4,10,-4:0.025:4,@(t)abs(sin(2*pi/4.*t)),'Seconds',false,'view(app.axh,[30,35])')
AbsSawtooth('正锯齿波',4,10,-4:0.025:4,@(t)abs(sawtooth(2*pi/4.*t,.745)),'Seconds',false,'view(app.axh,[30,35])')
Cos('余弦波',1,1,0:0.025:4,@(t)cos(2*pi.*t),'Seconds',false,'view(app.axh,[30,35])')
CosAddCos('余弦波相加波',1,11,0.5:0.002:2.5,@(t)(cos(2*pi*10.*t)+cos(2*pi*11.*t)),'Seconds',false,'view(app.axh,[30,35])')
Diric7('狄利克雷波_奇数',2*pi,7,linspace(-2*pi,2*pi,301),@(t)diric(t,7),'Rad',false,'view(app.axh,[60,70])')
Diric8('狄利克雷波_偶数',4*pi,7,linspace(-2*pi,2*pi,301),@(t)diric(t,8),'Rad',false,'view(app.axh,2)')
end
endclassdef FS_Visual3D_x_y < handle
%FS_VISUAL3D 可视化傅里叶级数，
% 用傅里叶复数系数重构信号函数，沿时间轴展开绘制任意周期函数3D曲线
% 绘制复数信号函数在实轴及虚轴上的投影
% 可视化动画加深对傅里叶级数的理解
%============================================================
% 一年多以前就想动手写这个代码，但没有想明白从哪入手
% 疫情期间不能上班正好可以静下来写代码
%============================================================
% 作者：易
% 964615453@qq.com
% 2022-4-23 8:36 完成2D的基本功能
% 2022-4-25 完成3D绘图功能
% 2022-4-26 加入实轴和虚轴的投影
% 2022-4-28 整合现有方法，增加最终接口完善本类
% 2022-4-29 修改bug，规范代码
% 2022-5-2 输入多种信号函数验证代码的正确性
% 2022-5-3 3:40 修正幅值为2倍，旋转方向，axLimt
properties
axh %图像对象句柄
Ln=[] %线对象数组，初始化为空
Rn=[] %圆对象数组，初始化为空，圆对象实际为线对象，方便3D绘图
RnCk=[] %保存每一个圆的圆心当前的坐标，复数
Lfun %线对象，最终绘制轨迹，傅里叶级数累加结果绘制的曲线
xLfun %线对象，实部投影
yLfun %线对象，虚部投影
fMsg Msg %广播帧数据消息对象
axLimt %xy轴范围限制
end
events
oneframe %广播一帧数据
end
%公共方法
methods
function obj = FS_Visual3D_x_y(axh)
if nargin == 1
obj.axh=axh;
end
obj.fMsg=Msg([]); %实例化事件消息对象
end
function FourierSeriesPlot(obj,C0,Cn,fn,t,SR,enRadio)
%C0 这个参数暂时不用，如果需要直流成分量可以加上，但旋转轴直线需要重新定位
%Cn 傅里叶级数系数向量数组
%fn 与傅里叶级数向量对应的频率数组，这个参数似乎有些多余，
% 但当初设计时考虑到Cn、fn是一一对应的，但可以不连续
%t 采样时间序列
%enRadio 是否广播图像帧，true、false
%SR 时间坐标，秒或弧度
init_ax(obj); %初始化轴对象
line(obj.axh,0*t,0*t,t); %画出时间线
axLimt = floor(sum(abs(Cn)).*2+.5)+1; %2倍幅值
xlim(obj.axh,[-axLimt,axLimt]);
ylim(obj.axh,[-axLimt,axLimt]);
obj.axLimt = axLimt;
switch SR
case 'Seconds'
%str=['Time' char(40) 'Seconds' char(41)];
zlabel(obj.axh,'Time(Seconds)');
case 'Rad'
%str=['Time' char(40) 'Rad' char(41)];
zlabel(obj.axh,'Time(Rad)');
end
N=length(t);
for k = 1:N
tk=t(k);
Ct=Cn.*exp(1j.*fn.*2.*pi.*tk).*2; %2倍幅值，tk时刻每谐波分量旋转向量的位置，复数
FS_Plot3D(obj,Ct,tk);
drawnow
if enRadio == true
obj.fMsg.Data=getframe(obj.axh); %获取axh轴对象的一帧数据
notify(obj,'oneframe',obj.fMsg); %广播一帧数据，选择存储Gif文件时触发
end
end
end
end
%私有方法
methods (Access = private)
function FS_Plot3D(obj,Cn,tk)
N = length(Cn); %傅里叶级数的个数
if isempty(obj.Ln) && isempty(obj.Rn)
%绘制第一帧并创建圆和线对象数组
Ci=0;
%这个初始值代表傅里叶系数中的C0项，即直流量，令其为零并不会改变波形
obj.RnCk = [obj.RnCk Ci];
for i= 1:N
[Ck,L_h,R_h] = RnLn_Plot3D(obj,i,Ci,Cn(i),tk);
Ci=Ck;
obj.Ln = [obj.Ln L_h]; %保存线对象
obj.Rn = [obj.Rn R_h]; %保存圆对象
obj.RnCk = [obj.RnCk Ck]; %保存当前圆心的坐标
end
x=[real(Ck) real(Ck)];
y=[imag(Ck) imag(Ck)];
z=[tk tk];
obj.Lfun = line(obj.axh,x,y,z, ...
'color','r','LineStyle','-','LineWidth',3); %绘制轨迹
axLimt = obj.axLimt;
xx=[-axLimt+0.1 -axLimt+0.1];
obj.xLfun = line(obj.axh,xx,y,z, ...
'color','k','LineStyle',':','LineWidth',1); %绘制实部投影线
yy=[axLimt-0.1 axLimt-0.1];
obj.yLfun = line(obj.axh,x,yy,z, ...
'color','b','LineStyle','-','LineWidth',1); %绘制虚部投影线
else
Ci=0;
for i= 1:N
Ck = RnLn_Update3D(obj,i,Ci,Cn(i),tk); %更新圆的位置及级数向量的位置
Ci=Ck;
end
x=get(obj.Lfun,'XData');
y=get(obj.Lfun,'YData');
z=get(obj.Lfun,'ZData');
x=[x real(Ck)];
y=[y imag(Ck)];
z=[z tk];
axLimt = obj.axLimt;
set(obj.Lfun,'XData',x,'YData',y,'ZData',z); %更新轨迹
xx=[get(obj.xLfun,'XData') -axLimt+0.1];
set(obj.xLfun,'XData',xx,'YData',y,'ZData',z); %更新实部投影线
yy=[get(obj.yLfun,'YData') axLimt-0.1];
set(obj.yLfun,'XData',x,'YData',yy,'ZData',z); %更新虚部投影线
end
end
function init_ax(obj)
%首次进入初始化轴对象
cla(obj.axh);
hold(obj.axh,"on");
axis(obj.axh,'equal');
box(obj.axh,'on');
xlabel(obj.axh,'Real');
ylabel(obj.axh,'Imag');
zlabel(obj.axh,'Time Axis');
title(obj.axh,'Fourie Series 3D Plot','FontSize',18,'FontWeight','bold');
end
function [Ck,L_h,R_h]= RnLn_Plot3D(obj,i,Ci,Cj,tK)
% 创建园对象和线对象，这里的圆也是线对象，使用line绘制的线对象
% i 第i次谐波
% Ci 给定相对坐标,复数,其为之前的1——i傅里叶级数之和
% Cj 第j傅里叶级数，复数
% Ck 输出下一点的坐标，复数，1——k傅里叶级数之和
% R_h 返回一个圆对象，便于用set移动圆的位置
% L_h 返回一个线对象，便于用set移动线的位置
% tk z轴坐标，时间,tk时刻
xi=real(Ci);
yi=imag(Ci);
xj=real(Cj);
yj=imag(Cj);
Ck=Ci+Cj;
L_h=line(obj.axh,[xi,xj+xi],[yi,yj+yi],[tK tK],'color','b','LineWidth',1); %绘制Ci向量线
rj=abs(Cj); %j次谐波所对应的圆半径，圆半径即为复数的摸
dCx=real(Ci); %圆心的x坐标
dCy=imag(Ci); %圆心的y坐标
w=0:0.04:2*pi;
if i == 1
x=rj.*cos(w); %1次谐波圆周各点，x数据，1次波不加偏移
y=rj.*sin(w); %1次谐波圆周各点，y数据，1次波不加偏移
else
x=rj.*cos(w)+dCx; %j次谐波圆周各点，x+偏移数据
y=rj.*sin(w)+dCy; %j次谐波圆周各点，y+偏移数据
end
z=ones(length(x),1).*tK; %时间坐标，j次谐波圆移动至tk
% 绘制园
R_h=line(obj.axh,x,y,z','LineStyle',':','LineWidth',0.2,'Color','b');
end
function Ck=RnLn_Update3D(obj,i,Ci,Cj,tk)
% 更新圆对象和线对象
xi=real(Ci);
yi=imag(Ci);
xj=real(Cj);
yj=imag(Cj);
Ck=Ci+Cj;
set(obj.Ln(i),'XData',[xi,xj+xi],'YData',[yi,yj+yi],'ZData',[tk tk]); %更新绘制Ci向量线
%rj=abs(Cj); %计算Cj向量所对应的圆半径
dC=Ci-obj.RnCk(i); %计算圆心偏移坐标
obj.RnCk(i)=Ci; %保存当前圆的圆心位置
dCx=real(dC);
dCy=imag(dC);
x=get(obj.Rn(i),'XData');
y=get(obj.Rn(i),'YData');
x=x+dCx;
y=y+dCy;
z=ones(length(x),1).*tk;
if Ci==0
set(obj.Rn(i),'ZData',z');
else
set(obj.Rn(i),'XData',x,'YData',y,'ZData',z'); %更新圆对象
end
end
end %methods
end %classdef
有需要完整代码的读友请将你的邮箱发给我，我的邮箱：964615453@qq.com