起源
之前制作钢琴频谱识别时候,发现傅里叶变换的神奇之处,激起了内心的好奇。前几天在bilibili又发现了用傅里叶变换绘图的视频。决定自己做一个,顺便仔细学习一下傅里叶变换。
过程
从对傅里叶变换一窍不通,到反复学习推导过程,以及自己边做边发现问题,加上自己的思考和参考视频中的讲解,终于对傅里叶变换有了一个比较深刻的理解。
原理
对于一个波,对其进行傅里叶变换之后可以得到若干角速度、初相角、大小不同的复数向量。
复数向量可以以箭头形式在复平面内表示。将这些箭头首尾相连,其运动一端的轨迹便可以在复平面内形成二维图形。
而此二维图形与输入的波之间有一定联系。
依据此关系,最终可以用一定的波进行傅里叶变换之后,使变换结果在复平面内绘出对应于原波的图形。
成果
演示动图:
程序:
ps:
1.需要注意的地方看程序注释。
2.输入图片需要是单连通的,具体定义可以百度。对于一些多连通的图形,用ps或画图 稍作修改即可成为单连通图形(比如演示中的“A” 和音乐符号)。
3.图形中尽量不要出现尖锐的角,由于寻找边界的算法很简陋(自己拍脑袋随便想了一下搞出来的),可能会寻到到尖锐处就无法继续。
function fft_draw()
%picture=imread('C:\Users\zxy\Pictures\Matlab\doo.png');
picture=imread('timg.png');
% 1 , 务必使用 白 底 黑 图 案 的 图片
% 2,因为只能绘制外轮廓,所以 务必使用 单 连 通 图形
% 3,图形务必居中( 大致居中即可 )。
N_point=800;%最终绘制时的采样点数,越大越精细,一般400足够
N_arrows=1;%显示的向量数;N_arrows=1或超过最大值时,自动取最大值
show=1;% show=1时显示寻找边界的动画,为其他值时不显示
cycle=2;% cycle=-1时无限循环显示动画,为n循环显示n次
picture=RuiHua(picture);%锐化
data_margin=ShiBie(picture, show);%识别边界
Show(data_margin,N_point, N_arrows,cycle);
function result=RuiHua(picture)%锐化,保证边界拾取的可靠性
[width, high,~]=size(picture);
for i=1:1:width
for j=1:1:high
if picture(i,j,1)<125
picture(i,j,:)=0;
else
picture(i,j,:)=255;
end
end
end
result=picture;
end
function data=ShiBie(picture, show)
%识别边界,得到边界点在图像中的坐标的数组(Xi,Yi),其中i=1,2,3...n
[width, high, ~]=size(picture);%获取图片宽度、高度
high_find=round(high/2);%从图片1/2高度开始寻找边界点
width_find=0;
%边界识别的起始点
for i=2:1:width-3
if abs( picture(high_find,i,1)- picture(high_find,i+1,1))>125 && width_find==0
%如果红色R(red)数组中,左侧点的亮度值 与 右侧点亮度值 的差 大于 125,则认为此处是边界。
%并记录此处坐标(high_find,width_find).记录的是边界上-内侧-的点的坐标
width_find=i+1;
end
end
margin=[high_find ,width_find];%边界点坐标数组,第一个点的坐标
picture_show=imshow(picture);
pause(2);
%开始识别边界
while picture(high_find,width_find ,3)-picture(high_find,width_find ,2)<100
%在经锐化的黑白图中,每个点RGB值只可能是(255,255,255)或(0,0,0).
%--如果出现例如(0,0,255)的情况,说明这个点已经被标记过。
%--又因为图形区域是单连通区域。再次遇到已经被标记过的点,说明所有边界已经识别完成
picture(high_find, width_find , 3)=255;
%将B(RGB的R)数组值改为(255),RG保持不变(0,0),作为已经被记录的边界点的标志。
if show==1
picture(high_find-1:high_find+1,width_find-1:width_find+1 , 1)=255;
%将识别点上下左右1的范围内的点的 R值全置为255,显示边界更明显
set(picture_show,'Cdata',picture);
pause(0.001);
end
for i=0:0.7854:5.5 % 0.7875 ~ pi/4
x_find_1=width_find+round(cos(i ) );
y_find_1=high_find+round( sin(i ) );
x_find_2=width_find+round(cos(i + 0.7854) );
y_fing_2=high_find+round( sin(i + 0.7854 ) );
%围绕上一次的识别点,从右侧开始,逆时针方向依次对比相邻两个点的 B值
%如果 B 值出现大幅下降,说明此处是边界,且较小值所在点是图形内侧。
if picture(y_find_1, x_find_1, 2) -picture(y_fing_2, x_find_2, 2) >125
high_find=y_fing_2;
width_find=x_find_2;
margin=[margin ;high_find,width_find];%记录新的边界点的坐标
end
end
end
data=margin;
end
function Show(data_margin, N_point, N_arrows,cycle)
[N,~]=size(data_margin);%边界点数
data_point=zeros(N_point,2);
arrows=plot([0 0],[0 0],'marker','.','color','b');hold on;%箭头
orbit=plot([0 0],[0 0],'r');%轨迹
xdata_X=zeros(1,2*N);
ydata_X=xdata_X;
xdata_Y=xdata_X;
ydata_Y=xdata_X;
data_margin=data_margin/100;
data_margin(: ,1)=data_margin(: ,1) - data_margin( 1,1);%平移图形,使起/止点坐标为(0,0)
data_margin(: ,2)=data_margin(: ,2) - data_margin( 1,2);
data_X=fft(data_margin(:,1));%对 x 傅里叶变换
data_Y=fft(data_margin(:,2));%对 y 傅里叶变换
[n, ~]=size(data_X);
ct=zeros(1,n);%各阶向量的初始角速度
limt=max(abs(data_X)*0.7);%控制绘出图形的显示大小
xlim([-limt limt]);
ylim([-limt limt]*0.8);
n=round(n/2);%阶数,傅里叶变换得到的是对称数组,取一半
if N_arrows==1 || N_arrows>n
N_arrows=n;
end
corner=zeros(n,2);
while cycle==-1 || cycle>0
if cycle>0
cycle=cycle-1;
end
for n_point=1:1:N_point
for i=1:1:n
ct(i)=ct(i)+(i-1)/N_point;%角速度
%理论上ct(i)= 0,pi, 2pi, 3pi ... npi
%在此比例为0:1:2:3:...n即可,具体大小只影响动画速度
if n_point==1 %--------------------------- x(t) ----------------------
if abs( data_X( i ) )>0
if imag(data_X( i) )>0
corner(i,1 )=acos( real(data_X(i))/abs(data_X(i)) );
%用傅里叶变换得到的复平面向量计算初始转角
else
corner( i ,1 )=-acos( real(data_X(i))/abs(data_X(i)) );
end
else
corner(i,1 )=0;
end
if abs( data_Y( i ) )>0 % ---------------------- y(t) ----------------------
if imag(data_Y( i) )>0%初始转角
corner(i ,2 )=acos( real(data_Y(i))/abs(data_Y(i)) );
else
corner(i ,2 )=-acos( real(data_Y(i))/abs(data_Y(i)) );
end
else
corner(i ,2 )=0;
end
end
if i>1% ------------- x(t) ---------------
%对于X而言,只需要x方向的运动,因此需要附加一个反向运动的向量 将y方向的运动抵消
xdata_X( i*2-1 )= 0.1*abs(data_X(i)) * cos( 2*pi*ct(i) + corner( i ,1 ) )+ xdata_X( i*2-2 );
xdata_X( i*2 )=0.1*abs(data_X(i)) * cos( -2*pi*ct(i) - corner( i ,1 ) )+ xdata_X( i*2-1);
ydata_X( i*2-1 )=0.1*abs(data_X(i)) * sin( 2*pi*ct(i) +corner( i ,1 ) ) + ydata_X( i*2 -2 );
ydata_X( i*2 )=0.1*abs(data_X(i)) * sin( -2*pi*ct(i) -corner( i ,1 ) ) + ydata_X( i*2-1);
else
xdata_X(i)=0.2*abs(data_Y(i)) * cos( corner( i ,1 ) );
ydata_X(i)=0;
end
if i>1% ------------y(t) ----------------
%对于Y而言,只需要y方向的运动,因此需要附加一个反向运动的向量 将x方向的运动抵消
xdata_Y( i*2-1 )=-0.1*abs(data_Y(i)) * sin( 2*pi*ct(i) + corner( i ,2 ) )+ xdata_Y( i*2-2 );
xdata_Y( i*2 )=-0.1*abs(data_Y(i)) * sin( -2*pi*ct(i) - corner( i ,2 ) )+ xdata_Y( i*2-1);
ydata_Y( i*2-1 )=0.1*abs(data_Y(i)) * cos( 2*pi*ct(i) +corner( i ,2 ) ) + ydata_Y( i*2-2 );
ydata_Y( i*2 )= 0.1*abs(data_Y(i)) * cos( -2*pi*ct(i) -corner( i ,2 ) ) + ydata_Y( i*2-1);
else
xdata_Y(i)=0;%
ydata_Y(i)=0.2*abs(data_Y(i)) * cos( corner( i ,2 ) );
end
end
x_draw=xdata_X(1:N_arrows) + xdata_Y(1:N_arrows) ;
y_draw=ydata_X(1:N_arrows)+ydata_Y(1:N_arrows);
data_point(n_point,:)=[y_draw(N_arrows) -x_draw(N_arrows)];
set(arrows,'xdata',y_draw(2:N_arrows) ,'ydata', -x_draw(2:N_arrows));
set(orbit,'xdata', data_point(1:n_point,1) ,'ydata', data_point(1:n_point,2) );
pause(0.009);
drawnow;
end
end
end
end