1701008931
1701008932
%最炫民族风 clc,clear,close all %清屏和清除变量 warning off %消除警告 fs = 44100; %采样率 dt = 1/fs; %采样时间 T16 = 0.125; t16 = [0:dt:T16]; [temp k] = size(t16); t4 = linspace(0,4*T16,4*k); %等分区间 t8 = linspace(0,2*T16,2*k); %等分区间 [temp i] = size(t4); %矩阵行、列提取 [temp j] = size(t8); %矩阵行、列提取 %Modification functions mod4=(t4.^4).*exp(-30*(t4.^0.5)); mod4=mod4*(1/max(mod4)); mod8=(t8.^4).*exp(-50*(t8.^0.5)); mod8=mod8*(1/max(mod8)); mod16=(t16.^4).*exp(-90*(t16.^0.5)); mod16=mod16*(1/max(mod16)); f0 = 2*146.8; %参考频率 ScaleTable = [2/3 3/4 5/6 15/16 … 1 9/8 5/4 4/3 3/2 5/3 9/5 15/8 … 2 9/4 5/2 8/3 3 10/3 15/4 4 … 1/2 9/16 5/8]; %1/4 notes do0f = mod4.*cos(2*pi*ScaleTable(21)*f0*t4); re0f = mod4.*cos(2*pi*ScaleTable(22)*f0*t4); mi0f = mod4.*cos(2*pi*ScaleTable(23)*f0*t4); fa0f = mod4.*cos(2*pi*ScaleTable(1)*f0*t4); so0f = mod4.*cos(2*pi*ScaleTable(2)*f0*t4); la0f = mod4.*cos(2*pi*ScaleTable(3)*f0*t4); ti0f = mod4.*cos(2*pi*ScaleTable(4)*f0*t4); do1f = mod4.*cos(2*pi*ScaleTable(5)*f0*t4); re1f = mod4.*cos(2*pi*ScaleTable(6)*f0*t4); mi1f = mod4.*cos(2*pi*ScaleTable(7)*f0*t4); fa1f = mod4.*cos(2*pi*ScaleTable(8)*f0*t4); so1f = mod4.*cos(2*pi*ScaleTable(9)*f0*t4); la1f = mod4.*cos(2*pi*ScaleTable(10)*f0*t4); tb1f = mod4.*cos(2*pi*ScaleTable(11)*f0*t4); ti1f = mod4.*cos(2*pi*ScaleTable(12)*f0*t4); do2f = mod4.*cos(2*pi*ScaleTable(13)*f0*t4); re2f = mod4.*cos(2*pi*ScaleTable(14)*f0*t4); mi2f = mod4.*cos(2*pi*ScaleTable(15)*f0*t4); fa2f = mod4.*cos(2*pi*ScaleTable(16)*f0*t4); so2f = mod4.*cos(2*pi*ScaleTable(17)*f0*t4); la2f = mod4.*cos(2*pi*ScaleTable(18)*f0*t4); ti2f = mod4.*cos(2*pi*ScaleTable(19)*f0*t4); do3f = mod4.*cos(2*pi*ScaleTable(20)*f0*t4); blkf = zeros(1,i); %1/8 notes do0e = mod8.*cos(2*pi*ScaleTable(21)*f0*t8); re0e = mod8.*cos(2*pi*ScaleTable(22)*f0*t8); mi0e = mod8.*cos(2*pi*ScaleTable(23)*f0*t8); fa0e = mod8.*cos(2*pi*ScaleTable(1)*f0*t8); so0e = mod8.*cos(2*pi*ScaleTable(2)*f0*t8); la0e = mod8.*cos(2*pi*ScaleTable(3)*f0*t8); ti0e = mod8.*cos(2*pi*ScaleTable(4)*f0*t8); do1e = mod8.*cos(2*pi*ScaleTable(5)*f0*t8); re1e = mod8.*cos(2*pi*ScaleTable(6)*f0*t8); mi1e = mod8.*cos(2*pi*ScaleTable(7)*f0*t8); fa1e = mod8.*cos(2*pi*ScaleTable(8)*f0*t8); so1e = mod8.*cos(2*pi*ScaleTable(9)*f0*t8); la1e = mod8.*cos(2*pi*ScaleTable(10)*f0*t8); tb1e = mod8.*cos(2*pi*ScaleTable(11)*f0*t8); ti1e = mod8.*cos(2*pi*ScaleTable(12)*f0*t8); do2e = mod8.*cos(2*pi*ScaleTable(13)*f0*t8); re2e = mod8.*cos(2*pi*ScaleTable(14)*f0*t8); mi2e = mod8.*cos(2*pi*ScaleTable(15)*f0*t8); fa2e = mod8.*cos(2*pi*ScaleTable(16)*f0*t8); so2e = mod8.*cos(2*pi*ScaleTable(17)*f0*t8); la2e = mod8.*cos(2*pi*ScaleTable(18)*f0*t8); ti2e = mod8.*cos(2*pi*ScaleTable(19)*f0*t8); do3e = mod8.*cos(2*pi*ScaleTable(20)*f0*t8); blke = zeros(1,j); %1/16 notes do0s = mod16.*cos(2*pi*ScaleTable(21)*f0*t16); re0s = mod16.*cos(2*pi*ScaleTable(22)*f0*t16); mi0s = mod16.*cos(2*pi*ScaleTable(23)*f0*t16); fa0s = mod16.*cos(2*pi*ScaleTable(1)*f0*t16); so0s = mod16.*cos(2*pi*ScaleTable(2)*f0*t16); la0s = mod16.*cos(2*pi*ScaleTable(3)*f0*t16); ti0s = mod16.*cos(2*pi*ScaleTable(4)*f0*t16); do1s = mod16.*cos(2*pi*ScaleTable(5)*f0*t16); re1s = mod16.*cos(2*pi*ScaleTable(6)*f0*t16); mi1s = mod16.*cos(2*pi*ScaleTable(7)*f0*t16); fa1s = mod16.*cos(2*pi*ScaleTable(8)*f0*t16); so1s = mod16.*cos(2*pi*ScaleTable(9)*f0*t16); la1s = mod16.*cos(2*pi*ScaleTable(10)*f0*t16); tb1s = mod16.*cos(2*pi*ScaleTable(11)*f0*t16); ti1s = mod16.*cos(2*pi*ScaleTable(12)*f0*t16); do2s = mod16.*cos(2*pi*ScaleTable(13)*f0*t16); re2s = mod16.*cos(2*pi*ScaleTable(14)*f0*t16); mi2s = mod16.*cos(2*pi*ScaleTable(15)*f0*t16); fa2s = mod16.*cos(2*pi*ScaleTable(16)*f0*t16); so2s = mod16.*cos(2*pi*ScaleTable(17)*f0*t16); la2s = mod16.*cos(2*pi*ScaleTable(18)*f0*t16); ti2s = mod16.*cos(2*pi*ScaleTable(19)*f0*t16); do3s = mod16.*cos(2*pi*ScaleTable(20)*f0*t16); blks = zeros(1,k); part0 = [mi1f la0e la0e do1f mi1f … re1e re1s mi1s re1e do1e re1e do1e la0f … mi1f la0e la0e do1f mi1f … so1e re1s mi1s re1e do1e re1e do1e ti0e so0e … mi1f la0e la0e do1f mi1f … re1e re1s mi1s re1e do1e re1e do1e la0e so0e … mi1f la0e la0e do1f mi1f … so1e mi1e blkf blkf blkf ]; part1 = [la0f la0e so0e la0f la0e do1e … do1f re1e do1e la0f la0f … do1f do1e so0e do1e re1e mi1e so1e … so1e mi1e re1f mi1f mi1f … la1e la1e la1e so1e mi1e mi1f do1e … la0e la0e la0e mi1e re1s mi1s re1e re1f … mi1e mi1e
1701008933
1701008934
运行程序，用户将听到《最炫民族风》歌曲。
1701008935
1701008936
音乐之所以和谐美妙，音乐就是一个数学的巧合，很大程度上得益于数学上的两个约等式同时成立：
1701008937
1701008938
2^（7/12）= 1.4983 ≈ 3/2，误差0.1％
1701008939
1701008940
2^（4/12）= 1.2599 ≈ 5/4，误差0.8％
1701008941
1701008942
1701008943
1701008944
1701008945
我和数学有约：趣味数学及算法解析 [:1701004240]
1701008946
我和数学有约：趣味数学及算法解析 7.4 玩转Google地球
1701008947
1701008948
谷歌地球（Google Earth，GE）是一款Google公司开发的虚拟地球仪软件，它把卫星照片、航空照相和GIS布置在一个地球的三维模型上。Google Earth于2005年向全球推出，被“PC世界杂志”评为2005年全球100种最佳新产品之一。用户们可以通过下载到自己电脑上的客户端软件，免费浏览全球各地的高清晰度卫星图片。
1701008949
1701008950
具体的Google地球如图7-15所示。
1701008951
1701008952
1701008953
1701008954
1701008955
图7-15 Google地球
1701008956
1701008957
手动旋转Google地球，得到中国板块视图，如图7-16所示。
1701008958
1701008959
1701008960
1701008961
1701008962
图7-16 Google地球之中国
1701008963
1701008964
采用程序进行驱动Google地球，设定视角，程序如下：
1701008965
1701008966
% 创建一个COM server，驱动Google Earth clc,clear,close all %清屏和清除变量 warning off %消除警告 handle = actxserver (‘googleearth.ApplicationGE’); %设定视角位置 lat = 37.77742; long = -122.4202; alt = 10000; altMode = 2; range = 0; tilt = 0; heading = 0; speed = 1; %设置Google Earth视角参数 handle.SetCameraParams(lat,long,alt,altMode,range,tilt,heading,speed); %获取视角中心目标点 my_target = handle.GetPointOnTerrainFromScreenCoords(0,0); pause(5) %设置聚焦点 new_target = handle.GetCamera(0); set(new_target, ‘FocusPointLatitude’, 37.8262); set(new_target, ‘FocusPointLongitude’, -122.4996); set(new_target, ‘FocusPointAltitude’, 235); set(new_target, ‘FocusPointAltitudeMode’, 2); set(new_target, ‘Range’, 0); set(new_target, ‘Tilt’, 75); set(new_target, ‘Azimuth’, 105); handle.SetCamera(new_target,1);
1701008967
1701008968
运行程序得到如图7-17所示的结果。
1701008969
1701008970
1701008971
1701008972
1701008973
图7-17 Google地球驱动
1701008974
1701008975
从Google地球中获取实际地形数据，可以用于地形地貌分析，这个操作显得尤为重要，借助于程序，用户需要了解地貌所在的位置，才能精确的进行截取，并得到地貌数据，具体的程序如下：
1701008976
1701008977
function ysw2_3 %怎么读取Google地球数据 %open GR ge_app = actxserver(‘GoogleEarth.ApplicationGE’); %GoogleEarth 6.0测试 %等待安装 pause(7) %设定视角 camPos = ge_app.GetCamera(0); %视角位置 camPos.FocusPointLatitude = 51.10719595658241; camPos.FocusPointLongitude = 13.638689585783112; camPos.FocusPointAltitude = 0.0; camPos.FocusPointAltitudeMode = ‘AbsoluteAltitudeGE’; camPos.Range = 11780.53549488341; camPos.Tilt = 0; camPos.Azimuth = 7; %0 < Speed <5 —> Speed> 5: Teleport ge_app.SetCameraParams(camPos.FocusPointLatitude, camPos.FocusPointLongitude, camPos.FocusPointAltitude, camPos. FocusPointAltitudeMode, camPos.Range, camPos.Tilt, camPos.Azimuth, 6); %measure x and y - distances u_li = ge_app.GetPointOnTerrainFromScreenCoords(-1,-1); %左下方 u_re = ge_app.GetPointOnTerrainFromScreenCoords(1,-1); %右下方 o_li = ge_app.GetPointOnTerrainFromScreenCoords(-1, 1); %左上方 dist_u = pos2dist(u_li.Latitude ,u_li.Longitude, u_re.Latitude , u_re. Longitude ,1); %distance x dist_o = pos2dist(u_li.Latitude ,u_li.Longitude, o_li.Latitude , o_li. Longitude ,1); %distance y %-1到1上产生网格 [X,Y] = meshgrid(-1:(1/30):1, -1:(1/30):1); Z = NaN(size(X)); %高度数据 %获取海拔 from GE for i= 1: 1 : size (X,1) for j =1:1: size (X,2) pot = ge_app.GetPointOnTerrainFromScreenCoords(X(i,j), Y(i,j) ); Z(i,j) = pot.Altitude; end; end; %计算实际高度 X_dist = (X + 1).*0.5.*dist_u; %km Y_dist = (Y + 1).*0.5.*dist_o; %km Z_dist = Z; %海拔m %获取RGB图像 from GE jmouseemu([600,500],‘normal’); %设置视角 inputemu(‘key_alt’,‘c’); %复制图像 pause(2); %等待 [imdata] = clipboardimage(‘jpg’); %得到图像 %转换为索引图像 [X2, map2] = rgb2ind(imdata, 255); %uint8 %绘制截取曲面 figure(‘Color’, [1 1 1 ]); surface(X_dist,Y_dist,Z_dist,flipud(X2),‘FaceColor’,‘texturemap’,‘EdgeColor’,‘none’, ‘CDataMapping’,‘direct’) colormap(map2); xlim([0 dist_u]); ylim([0 dist_o]); view([-12,72]); box on grid on
1701008978
1701008979
运行程序得到如图7-18所示的结果。
1701008980