傅里叶变换的积分实现
过冷水最近遇到了这么一个问题。如何正确实现上图所表示的图像函数相互转换。可以看出图像图像很复杂，用一般的函数并不能准确的去描述图像。至于图像的转换公式，天！复杂的积分公式，理论描述该问题是如此的简单，过冷水往期也有和大家一起分享复杂函数的积分问题，本期过冷水会带大家一起做一下两幅图像的相互转换工作，重点讲一下积分计算中的小技巧。
起初冷水拿到手的数据是这样的：
%数据
x1=[2.60252365930599,2.67305657466650,2.74344115582014,2.76635879077955,2.81308406593950,2.83555669827833,2.85713932537603,2.87827694985315,2.90052708088170,2.92210970797939,2.94391483638737,2.96683247134678,2.98804426292733,3.00947855581816,3.03158035263985,3.02950367374386,3.07526477655924,3.16908616239629,3.24006408037736,3.26446505740534,3.26557756395677,3.26631923499105,3.29138771594988,3.29346439484588,3.31793953897729,3.36592565489552,3.39099413585436,3.41546927998576,3.46397456562800,3.51203484864966,3.53658415988450,3.58494111131988,3.63314972854840,3.72919612748831,3.82502002511793,3.94420656032752,4.01585198223945,4.13474184903533,4.30072782650831,4.39558755179336,4.53769172196236,4.72666950149818,4.91535061262028,5.03342464127844,5.12776519683949,5.22188325109026,5.33958644423128,5.41019352669521,5.50431158094598,5.62223727539729,5.71628116254462,5.81047338389882,5.88108046636275,6.02273963391118,6.16462130276989,6.30679964004233,6.42539283842449,6.54420853811695,6.66317257201626,6.80572174480583,6.92438911029143,7.04298230867360,7.16150133995234,7.30367967722477,7.42205037429665,7.61132482224618,7.77664329578830,7.96569524242754,8.13108788307309,8.24901357752440,8.43784302285336,8.67428774858341,8.88707316831977,9.12351789404982,9.36040762240044,9.57348971055052,9.73947568802350,9.95270611038043,10.1657881985305,10.4734333435519,10.7808559872630,10.9699079339023,11.1590340476449,11.3954787733750,11.6083383602148,11.8685165590420,12.1761617040634,12.3891696251100,12.6022517132601,12.8625040791907,13.1938085302058,13.5487722872146,13.7617060411578,14.1639142430505,14.3531886910001,14.8026413377768];
x2=[0.979407169516535,1.11285912681191,1.31297196423922,1.44622862598733,1.57948528773544,1.66788907212916,1.75590226542835,1.77764516969013,1.79919277840465,1.84306917802274,1.86403090009548,1.88544831177849,1.90699592049301,1.92860862772329,1.94937505424876,1.97053207186876,1.99168908948876,2.01336689523479,2.03387292769725,2.05522524086451,2.07638225848451,2.09806006423054,2.14154587275410,2.23001475566357,2.27499783004948,2.29830309868935,2.34387205971704,2.36698203280965,2.39048259699679,2.41411335821543,2.45968231924312,2.48292248936724,2.50583716691259,2.52888204148945,2.57438590400139,2.59730058154674,2.66500303793074,2.79930127593091,2.95592830483465,3.06744206232098,3.20115441367937,3.35673986633105,3.48993142956341,3.57833521395712,3.68887249370714,3.73255359777797,3.86535456991581,3.99809044353789,4.06455602812256,4.24253537019356,4.30952174290426,4.42136099296936,4.51093655064665,4.69015276451697,4.82432080548564,4.98042704626335,5.13646818852530,5.40337210311605,5.60328964499609,5.73641610971270,5.98073083933686,6.13625119347279,6.24756965541186,6.38141220380175,6.55984723548303,6.69381998090444,6.73847756271157,7.02810085930041,7.18407690304661,7.40658362989324,7.60682666435205,7.76234701848798,7.94026126204323,8.09597691172641,8.27408645082892,8.38547001128374,8.56370974741776,8.74194948355178,8.96471660446142,9.16528513149900,9.36559326447357,9.54357260654457,9.74381564100339,9.94418887249371,10.1666955993403,10.3448702369586,10.5675071608367,10.7679454908428,10.9683838208489,11.2133495356306,11.3913939762173,11.6139658015797,11.7920102421665,12.0145169690131];
y1=[0,0.0564263322884013,0.131661442006270,0.225705329153605,0.300940438871473,0.451410658307210,0.714733542319749,1.03448275862069,1.21316614420063,1.47648902821317,1.71159874608150,1.80564263322884,2.11598746081505,2.39811912225705,2.59561128526646,2.85893416927900,3.05642633228840,3.15987460815047,3.15987460815047,3.06583072100313,2.92476489028213,2.83072100313480,2.65203761755486,2.38871473354232,2.28526645768025,2.20062695924765,2.02194357366771,1.91849529780564,1.76802507836991,1.67398119122257,1.56112852664577,1.42946708463950,1.31661442006270,1.13793103448276,0.987460815047022,0.874608150470219,0.789968652037618,0.714733542319749,0.667711598746082,0.639498432601881,0.620689655172414,0.658307210031348,0.733542319749216,0.761755485893417,0.799373040752351,0.865203761755486,0.940438871473354,0.987460815047022,1.05329153605016,1.10031347962382,1.17554858934169,1.23197492163009,1.27899686520376,1.31661442006270,1.32601880877743,1.29780564263323,1.26018808777429,1.19435736677116,1.10971786833856,1.03448275862069,0.987460815047022,0.949843260188088,0.921630094043887,0.893416927899687,0.884012539184953,0.884012539184953,0.921630094043887,0.949843260188088,0.978056426332288,1.02507836990596,1.08150470219436,1.10031347962382,1.11912225705329,1.13793103448276,1.10031347962382,1.08150470219436,1.03448275862069,0.996865203761756,0.978056426332288,0.968652037617555,0.987460815047022,1.01567398119122,1.03448275862069,1.05329153605016,1.06269592476489,1.07210031347962,1.06269592476489,1.05329153605016,1.03448275862069,1.03448275862069,1.02507836990596,1.01567398119122,1.01567398119122,1.01567398119122,1.01567398119122,1.02507836990596];
y2=[0.0248839076469057,0.0422326620953042,0.0770278187657321,0.120692648207621,0.164357477649510,0.251947530596302,0.392169733530076,0.462313384254839,0.558773109973093,0.646428261435639,0.821836212134363,0.935839987848277,1.03229971356653,1.11998741428695,1.32171143997917,1.47080331568440,1.61989519138964,1.69881086711223,1.93562299279576,2.05839879350751,2.20749066921274,2.28640634493534,2.42669364638486,2.50551167433383,2.44404240083326,2.30365745161010,2.16323995312907,2.04917107889940,1.88247005468275,1.69822498047045,1.55780748198941,1.42619455776408,1.33844175852791,1.23314490929607,1.10149943581286,1.01374663657669,0.890840638833435,0.794153068310043,0.688660923530943,0.662182102248069,0.644442756705147,0.679303011891329,0.731739866331048,0.819329919277841,0.924431472962417,1.03840269941845,1.14347170384515,1.25731273326968,1.30107521048520,1.31835886641785,1.29194514365072,1.22160619737870,1.15129980036455,1.00191498133843,0.922771460810694,0.887455516014235,0.860911596215606,0.895609105112404,0.956720336776322,1.01792921621387,1.09651939935769,1.14015167954171,1.13998893325232,1.10470553771374,1.06058501866157,1.00775757312733,0.990148424615919,0.963409209270029,0.945637314469230,0.962855871886121,0.980106978560889,1.02373925874490,1.04979493967538,1.06711114486590,1.06685075080288,1.05791597951567,1.04011153545699,1.02230709139832,1.00443754882389,0.977828530509504,0.986307612186442,1.00359126811909,1.02084237479385,1.02054943147296,1.03776798888985,1.02873556982901,1.02841007725024,1.01934510893152,1.01028014061279,1.00115007377832,1.00966170471313,1.01810823713219,1.02661986806701,1.04383842548390];用这样一组数据做积分运算不是太合适，我们首先的得到能够准确表述这些点的函数解析式。经过过冷水的反复测试，找到了有这几种类型函数可描述实点：'smoothingspline'、'gauss'
'lnterpolant'。
由图像可知上述几种函数都可以准确的得到函数图像，三种拟合函数各有其特点，gauss8可以得到明确的解析式，'smoothingspline'、'lnterpolant'也得到了拟合函数，奈何过冷水水平有限不知该如何获得可以看的解析式而不是打包函数。根据其函数类型不同自然求积分的方法也就不同了。过冷水先来讲一讲如何利用 'smoothingspline'、'lnterpolant'进行积分和函数转换。由于'smoothingspline'、'lnterpolant'函数不可以查看解析式，所以只能用大数定理或者蒙特卡洛来求积分:
由图像可发现用'smoothingspline'、'lnterpolant'做拟合再用大数定理求积分基本上能够得到我们想要转换效果看上去转换函数和原始值有 误差是由实验数据的不完全匹配引起的，我们使用的方法是没有问题的。
我们再用解析式的方式求积分我们用gauss8得到的两个解析式分别是：
g(y) | f(x) | |||||
|---|---|---|---|---|---|---|
ag | bg | cg | af | bf | cf | |
1 | 2.657 | 3.131 | 0.221 | 2.101 | 2.184 | 0.287 |
2 | 1302 | 3.532 | 0.472 | 1.121 | 4.121 | 0.749 |
3 | 1.228 | 6.102 | 1.717 | 0.879 | 6.064 | 1.232 |
4 | 0.709 | 8.827 | 1112 | 0.473 | 8.033 | 0.945 |
5 | 0.678 | 1.861 | 1.121 | 1.725 | 18.281 | 8.566 |
6 | 0.704 | 11.011 | 2.175 | 0.6928 | 2.648 | 0.605 |
7 | 0.9951 | 15.112 | 3.713 | 0.358 | 9.712 | 1.521 |
8 | 0.254 | 7.886 | 0.715 | 0.83 | 1.559 | 0.255 |
我们已经知道符号积分函数int,根据最原始的程序我们可以得到转换函数：
符号运算的代码：
g2= fit( x1, y1,'gauss8');f2= fit( x2, y2,'gauss8');
%g(y)转换为f(x)
f_x=linspace(min(x2),max(x2),50);
for i=1:length(f_x)
syms x
xgauss8=g2.a1.*exp(-((x-g2.b1)./g2.c1).^2) + g2.a2.*exp(-((x-g2.b2)./g2.c2).^2) + g2.a3.*exp(-((x-g2.b3)./g2.c3).^2) + g2.a4.*exp(-((x-g2.b4)./g2.c4).^2) + g2.a5.*exp(-((x-g2.b5)./g2.c5).^2) + g2.a6.*exp(-((x-g2.b6)./g2.c6).^2) + g2.a7.*exp(-((x-g2.b7)./g2.c7).^2) + g2.a8.*exp(-((x-g2.b8)./g2.c8).^2);
%xgauss8=fdis2.a1.*exp(-((x-fdis2.b1)./fdis2.c1).^2) + fdis2.a2.*exp(-((x-fdis2.b2)./fdis2.c2).^2) + fdis2.a3.*exp(-((x-fdis2.b3)./fdis2.c3).^2) + fdis2.a4.*exp(-((x-fdis2.b4)./fdis2.c4).^2) + fdis2.a5.*exp(-((x-fdis2.b5)./fdis2.c5).^2) + fdis2.a6.*exp(-((x-fdis2.b6)./fdis2.c6).^2) + fdis2.a7.*exp(-((x-fdis2.b7)./fdis2.c7).^2) + fdis2.a8.*exp(-((x-fdis2.b8)./fdis2.c8).^2);
int_xgauss8=@(x) sin(x.*f_x(i)).*x.*(xgauss8-1);
f_gauss8(i)=vpa((1+(4*pi)./(25*f_x(i))*int(int_xgauss8,x,min(x1),max(x1))),4);
end
%f(x)转换为g(y)
g_y=linspace(min(x1),max(x1),50);
for i=1:length(str_y)
syms x
ygauss8=f2.a1.*exp(-((x-f2.b1)./f2.c1).^2) + f2.a2.*exp(-((x-f2.b2)./f2.c2).^2) + f2.a3.*exp(-((x-f2.b3)./f2.c3).^2) + f2.a4.*exp(-((x-f2.b4)./f2.c4).^2) + f2.a5.*exp(-((x-f2.b5)./f2.c5).^2) + f2.a6.*exp(-((x-f2.b6)./f2.c6).^2) + f2.a7.*exp(-((x-f2.b7)./f2.c7).^2)+ f2.a8.*exp(-((x-f2.b8)./f2.c8).^2);
int_ygauss8=@(x) sin(x.*g_y(i)).*x.*(ygauss8-1);
g_gauss8(i)=vpa((1+(1*35)./(2*pi.^2*g_y(i))*int(int_ygauss8,x,min(x2),max(x2))),4);
end观察上述代码，我们在计算的时候有使用for循环，之前过冷水一直觉得该语句OK!没有问题，随着学习的升入就会有想去有优化它的欲望，毕竟简洁是美！本案例可以用arrayfun函数代替循环于是就得到下段代码：
sym x
str_y=linspace(min(x1),max(x1),20);dis_x=linspace(min(x2),max(x2),2);
ygauss8=fstr2.a1.*exp(-((x-fstr2.b1)./fstr2.c1).^2) + fstr2.a2.*exp(-((x-fstr2.b2)./fstr2.c2).^2) + fstr2.a3.*exp(-((x-fstr2.b3)./fstr2.c3).^2) + fstr2.a4.*exp(-((x-fstr2.b4)./fstr2.c4).^2) + fstr2.a5.*exp(-((x-fstr2.b5)./fstr2.c5).^2) + fstr2.a6.*exp(-((x-fstr2.b6)./fstr2.c6).^2) + fstr2.a7.*exp(-((x-fstr2.b7)./fstr2.c7).^2)+ fstr2.a8.*exp(-((x-fstr2.b8)./fstr2.c8).^2);
F=@(x) sin(x.*str_y).*x.*(ygauss8-1)
dis_gauss8=vpa(arrayfun(@(str_y) eval(char(1+(1*35)/(2*pi.^2*str_y)*int(F,x,min(x2),max(x2)))),str_y,'UniformOutput',false),4);
xgauss8=fdis2.a1.*exp(-((x-fdis2.b1)./fdis2.c1).^2) + fdis2.a2.*exp(-((x-fdis2.b2)./fdis2.c2).^2) + fdis2.a3.*exp(-((x-fdis2.b3)./fdis2.c3).^2) + fdis2.a4.*exp(-((x-fdis2.b4)./fdis2.c4).^2) + fdis2.a5.*exp(-((x-fdis2.b5)./fdis2.c5).^2) + fdis2.a6.*exp(-((x-fdis2.b6)./fdis2.c6).^2) + fdis2.a7.*exp(-((x-fdis2.b7)./fdis2.c7).^2) + fdis2.a8.*exp(-((x-fdis2.b8)./fdis2.c8).^2);
G=@(x) sin(x.*dis_x).*x.*(xgauss8-1)
str_gauss8=vpa(arrayfun(@(dis_x) eval(char(1+(4*pi)/(25*dis_x)*int(F,x,min(x2),max(x2)))),str_y,'UniformOutput',false),4);基本上代码已经完善好了，but不完美之处在于符号求积分的效率太低，程序运行时间比较长，过冷水能和大家分享程序就不会把不完美的程序和大家分享出来一定会有好的解决办法。用quad函数可以迅速提高运算效率，but quad函数也是不完美的，是在积分的时候有限制。感兴趣的可以摸索一下问题在什么地方。
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2021-02-15，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
腾讯云开发者
