《2022年2022年离散信号MATLAB频谱分析程序 .pdf》由会员分享，可在线阅读，更多相关《2022年2022年离散信号MATLAB频谱分析程序 .pdf（14页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、离散信号MATLAB频谱分析程序【ZZ】%FFT 变换,获得采样数据基本信息，时域图，频域图%这里的向量都用行向量，假设被测变量是速度，单位为m/s clear; close all; load data.txt %通过仪器测量的原始数据，存储为data.txt 中，附件中有一个模版(该信号极不规则 ) A=data; %将测量数据赋给A，此时 A 为 N 2 的数组x=A(:,1); %将 A 中的第一列赋值给x，形成时间序列x=x; %将列向量变成行向量y=A(:,2); %将 A 中的第二列赋值给y，形成被测量序列y=y; %将列向量变成行向量%显示数据基本信息fprintf(n 数据基
2、本信息： n) fprintf( 采样点数= %7.0f n,length(x) %输出采样数据个数fprintf( 采样时间= %7.3f sn,max(x)-min(x) %输出采样耗时fprintf( 采样频率= %7.1f Hzn,length(x)/(max(x)-min(x) %输出采样频率fprintf( 最小速度= %7.3f m/sn,min(y) %输出本次采样被测量最小值fprintf( 平均速度= %7.3f m/sn,mean(y) %输出本次采样被测量平均值fprintf( 速度中值= %7.3f m/sn,median(y) %输出本次采样被测量中值fprintf
3、( 最大速度= %7.3f m/sn,max(y) %输出本次采样被测量最大值fprintf( 标准方差= %7.3f n,std(y) %输出本次采样数据标准差fprintf( 协 方 差 = %7.3f n,cov(y) %输出本次采样数据协方差fprintf( 自相关系数= %7.3f nn,corrcoef(y) %输出本次采样数据自相关系数%显示原始数据曲线图（时域）subplot(2,1,1); plot(x,y) %显示原始数据曲线图axis(min(x) max(x) 1.1*floor(min(y) 1.1*ceil(max(y) %优化坐标，可有可无xlabel( 时间 (
4、s);ylabel( 被测变量 y); title( 原始信号 (时域); grid on; %傅立叶变换y=y-mean(y); %消去直流分量，使频谱更能体现有效信息Fs=2000; %得到原始数据data.txt 时，仪器的采样频率。其实就是length(x)/(max(x)-min(x); N=10000; %data.txt 中的被测量个数，即采样个数。其实就是length(y); z=fft(y); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 14 页 -
5、 - - - - - - - - %频谱分析f=(0:N-1)*Fs/N;Mag=2*abs(z)/N; %幅值，单位同被测变量y Pyy=Mag.2; %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).2=X.2，故这里有很多表达方式%显示频谱图 (频域) subplot(2,1,2) plot(f(1:N/2),Pyy(1:N/2),r) %显示频谱图% | % 将这里的 Pyy 改成 Mag 就是 幅值-频率图了axis(min(f(1:N/2) max(f(1:N/2) 1.1*floor(min(Pyy(1:N/2) 1.1*ceil(max(Pyy(1:N/2
6、) xlabel( 频率 (Hz) ylabel( 能量) title( 频谱图 (频域) grid on; %返回最大能量对应的频率和周期值a b=max(Pyy(1:N/2); fprintf(n 傅立叶变换结果： n) fprintf( FFT_f = %1.3f Hzn,f(b) % 输出最大值对应的频率fprintf( FFT_T = %1.3f sn,1/f(b) %输出最大值对应的周期*% FFT实践及频谱分析 % %*% %*%*1.正弦波 *% fs=100;%设定采样频率N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率%生成正弦信号x=sin(
7、2*pi*f0*t); figure(1); subplot(231); plot(t,x);%作正弦信号的时域波形xlabel(t); ylabel(y); title(正弦信号y=2*pi*10t时域波形 ); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 14 页 - - - - - - - - - grid; %进行 FFT 变换并做频谱图y=fft(x,N);%进行 fft变换mag=abs(y);%求幅值f=(0:length(y)-1)*fs/length
8、(y);%进行对应的频率转换figure(1); subplot(232); plot(f,mag);%做频谱图axis(0,100,0,80); xlabel(频率 (Hz); ylabel(幅值 ); title(正弦信号y=2*pi*10t幅频谱图 N=128); grid; %求均方根谱sq=abs(y); figure(1); subplot(233); plot(f,sq); xlabel(频率 (Hz); ylabel(均方根谱 ); title(正弦信号y=2*pi*10t均方根谱 ); grid; %求功率谱power=sq.2; figure(1); subplot(234
9、); plot(f,power); xlabel(频率 (Hz); ylabel(功率谱 ); title(正弦信号y=2*pi*10t功率谱 ); grid; %求对数谱ln=log(sq); figure(1); subplot(235); plot(f,ln); xlabel(频率 (Hz); ylabel(对数谱 ); title(正弦信号y=2*pi*10t对数谱 ); grid; %用 IFFT恢复原始信号xifft=ifft(y); magx=real(xifft); ti=0:length(xifft)-1/fs; 名师资料总结 - - -精品资料欢迎下载 - - - - -
10、- - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 14 页 - - - - - - - - - figure(1); subplot(236); plot(ti,magx); xlabel(t); ylabel(y); title(通过 IFFT 转换的正弦信号波形); grid; %*2.矩形波 *% fs=10;%设定采样频率t=-5:0.1:5; x=rectpuls(t,2); x=x(1:99); figure(2); subplot(231); plot(t(1:99),x);%作矩形波的时域波形xlabel(t); yl
11、abel(y); title(矩形波时域波形); grid; %进行 FFT 变换并做频谱图y=fft(x);%进行 fft变换mag=abs(y);%求幅值f=(0:length(y)-1)*fs/length(y);%进行对应的频率转换figure(2); subplot(232); plot(f,mag);%做频谱图xlabel(频率 (Hz); ylabel(幅值 ); title(矩形波幅频谱图); grid; %求均方根谱sq=abs(y); figure(2); subplot(233); plot(f,sq); xlabel(频率 (Hz); ylabel(均方根谱 ); ti
12、tle(矩形波均方根谱); grid; %求功率谱power=sq.2; figure(2); subplot(234); plot(f,power); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 14 页 - - - - - - - - - xlabel(频率 (Hz); ylabel(功率谱 ); title(矩形波功率谱 ); grid; %求对数谱ln=log(sq); figure(2); subplot(235); plot(f,ln); xlabel(频
13、率 (Hz); ylabel(对数谱 ); title(矩形波对数谱 ); grid; %用 IFFT恢复原始信号xifft=ifft(y); magx=real(xifft); ti=0:length(xifft)-1/fs; figure(2); subplot(236); plot(ti,magx); xlabel(t); ylabel(y); title(通过 IFFT 转换的矩形波波形); grid; %*3.白噪声 *% fs=10;%设定采样频率t=-5:0.1:5; x=zeros(1,100); x(50)=100000; figure(3); subplot(231); p
14、lot(t(1:100),x);%作白噪声的时域波形xlabel(t); ylabel(y); title(白噪声时域波形); grid; %进行 FFT 变换并做频谱图y=fft(x);%进行 fft变换mag=abs(y);%求幅值f=(0:length(y)-1)*fs/length(y);%进行对应的频率转换figure(3); subplot(232); plot(f,mag);%做频谱图xlabel(频率 (Hz); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页
15、，共 14 页 - - - - - - - - - ylabel(幅值 ); title(白噪声幅频谱图); grid; %求均方根谱sq=abs(y); figure(3); subplot(233); plot(f,sq); xlabel(频率 (Hz); ylabel(均方根谱 ); title(白噪声均方根谱); grid; %求功率谱power=sq.2; figure(3); subplot(234); plot(f,power); xlabel(频率 (Hz); ylabel(功率谱 ); title(白噪声功率谱 ); grid; %求对数谱ln=log(sq); figure
16、(3); subplot(235); plot(f,ln); xlabel(频率 (Hz); ylabel(对数谱 ); title(白噪声对数谱 ); grid; %用 IFFT恢复原始信号xifft=ifft(y); magx=real(xifft); ti=0:length(xifft)-1/fs; figure(3); subplot(236); plot(ti,magx); xlabel(t); ylabel(y); title(通过 IFFT 转换的白噪声波形); grid 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - -
17、 - 名师精心整理 - - - - - - - 第 6 页，共 14 页 - - - - - - - - - 一个计算图像分形维数的软件这个是计算图象的分形维数的软件,可在 WINDOWS 下安装运行 ,运行时 ,新建一个project, 在其中导入图象即可下载地址：http:/ 程序http:/ http:/看下里面的留言区，有些东西哦一个求盒维数的Matlab 程序：%本程序只能计算时间序列的盒维数%box11 clear data1; % 存放时间序列的数据文件，共6 列c=0; n=210; ll=1:1383-1220+1; b(:,1)=ll; b(:,2)=aaa(1220:13
18、83,6)*100; %aaa是 data1 中的函数，里面是6 列时间序列%划分不同大小的格子，并统计不同量级格子中含有数据的格子数k=2; for p=1:9 .k=2p .stepx=(max(b(:,1)-min(b(:,1)+1)/k; .stepy=(max(b(:,2)-min(b(:,2)+1)/k; .mm(p)=0; .for i=1:k .i .for j=1:k .lowx=stepx*(i-1)+min(b(:,1)-0.5; .upx =stepx*i+min(b(:,1)-0.5; .lowy=stepy*(j-1)+min(b(:,2)-0.5; .upy =s
19、tepy*j+min(b(:,2)-0.5; .aa=find(b(:,1)=lowx & b(:,1)=lowy & b(:,2)upy); .if isempty(aa) .mm(p)=mm(p)+1; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 14 页 - - - - - - - - - .end .end .end end %在双对数坐标系中画点和曲线x(1)=1; for i=2:9 x(i)=2*x(i-1); end x=log(x); y=log(m
20、m); plot(x,y,*) hold on; v=polyfit(x(1:6),y(1:6),1); yy=polyval(v,x); plot(x,yy); hold off; title(v(1); 关联维数程序：% estimation of the correlation dimension % Grassberger and Procaccia algorithm r=0.50; m=22; t=5; dr=0.10; keyboard; % input r,m,t,dr n=length(fname); mm=n-(m-1)*t; eeaa=zeros(mm,m); for i
21、=1:mm for j=1:m eeaa(i,j)=fname(i+(j-1)*t); end; end; fenmu=mm*(mm-1); mmm=250/r; eebb=zeros(mmm,4); h=1; for ij=r:r:250 k=0; 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页，共 14 页 - - - - - - - - - for i=1:mm-1 for j=i+1:mm hh=norm(eeaa(i,-eeaa(j,inf); %(eeaa(i,:
22、）-eeaa(j,: ）,inf); if hhij k=k+2; end; end; end; crm=k/fenmu; eebb(h,1)=k; eebb(h,2)=crm; eebb(h,3)=log2(crm); eebb(h,4)=log2(ij); h=h+1; end; 建议用这个程序的人还是弄懂关联维数的计算理论！其实你弄懂了，不用我给你程序，你自己也可以编写出来！我在网上下载了一段求广义分形维数的MATLAB程序，如下：function dq,rq=fdim(q,x,trace); % from fdim.sh % Estimates generalised fractal
23、dimensions of a set of points x;% e.g. q=0 is Hausdorf, q=1 is Information, q=2 is Correlation, etc. % v0.20 1/7 /94, copyright AJ Roberts, arobertsusq.edu.au % Made publically available under the GNU arrangements.% % x(j,:) = coordinates of the jth point of a set; % should be embeded; % q = row vec
24、tor of exponents % dq = corresponding row of qth generalised dimension % rq = optimistic range-factor over which the dim holds % trace = optional argument: set to 1 if you want trace prints and graphics % otherwise omit or set to 0 % % The following parameters may be changed-at your own risk. rrat=2
25、; % square of radius-ratio for covering disks. tol=0.03; % approx error in minimisation of curve fit. if nargin0)/lrat ); nor(j,max(k)=0; for l=k; nor(j,l)=nor(j,l)+1; end;if unwarn, et=etime(clock,t0); if (et5)&(2*j et=(n-j)*5/j;disp(num2str(et) secs to go);unwarn=0;end; end; end; nor=cumsum(ones(1
26、,n);flipud(nor)/n; j,nr=size(nor); lrat=lrat/2;lrad=log(rmax)+(-nr+1):0)*lrat; % % fit the curves for each exponent q % Ref: Pawelzik & Schuster , Phys Rev A 35 (1987) pp481-484.lrad=lrad(1)-lrat lrad lrad(nr)+lrat; % pad out dq=; rq=; for qi=q, % if trace, disp(analysing q = num2str(qi); end;if abs
27、(qi-1)sqrt(eps) , lc=log( sum( nor.(qi-1) )/n )/(qi-1);else lc=sum( log(nor) )/n; end; j=find(lc0.9*lc(1);p=lrad(min(j) lrad(max(j)+2) 0.5 0.5; % p=fmins(rampdiff,p,0 tol,lrad,lc(1) lc 0,trace);p=fminsearch(rampdiff,p,lrad,lc(1) lc 0,trace); di=(lc(nr)-lc(1)/(p(2)-p(1);dq=dq di;rq=rq p(2)-p(1);if tr
28、ace,title(q = num2str(qi) Dq = num2str(di); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 14 页 - - - - - - - - - % pause(1); end % Dq(i)=dq; end; rq=exp(rq); % end;% END_OF_FILE % if test 2107 -ne wc -c hint.m feval Undefined function rampdiff.Error in = D:ma
29、tlab6.5toolboxmatlabfunfunfminsearch.m On line 125 = fv(:,1) = feval(funfcn,x,varargin:); Error in = D:51622470falphafdim.m On line 63 = p=fminsearch(rampdiff,p,lrad,lc(1) lc 0,trace);0 我要推荐Matlab 在分形模拟上的一些应用2008-10-21 09:21:25| 分类：默认分类| 标签：无 |字号大中小 订阅把在形态、功能和信息等方面具有自相似性的对象称为分形。自相似性是指局部的形态与整体的形态相似，局
30、部与整体相互依赖。由此可以通过局部认识和推及整体。在自然界中，很多的自然景观就具有自相似性。如云彩、山脉、海岸线、火焰、名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页，共 14 页 - - - - - - - - - 水波等，只要抽象出这些自然景观的某些特征，再不断放大，就可以得到整体。计算机分形图形就是利用这一原理实现的。MATLAB 作为一种科学计算软件，具有输入简捷、运算高效、内容丰富、扩展性强等特点，并且它还可与Word 结合起来，使含有 MATLAB 运行内容的文
31、章的编辑变得非常方便。这些特点使MATLAB 深受科技人员，特别是非软件专业的人们的喜爱。分形在所有尺寸下都有无限的细节，因而完全由计算机不可能产生精确的分形，只能给出近似的描绘，即模拟。而这里给出的程序中均有一个控制精度的量，只要计算机运算速度和屏幕分辨率允许，精度可任意提高。分形可作为双曲迭代函数系统的吸引子，根据程序中迭代的选取，将分形模拟分为确定迭代法和随机迭代法。Cantor 集选取一个欧氏长度的直线段，将该线段三等分，去掉中间一段，剩下两段。将剩下的两段分别再三等分，各去掉中间一段，剩下四段。将这样的操作继续下去，直到无穷，则可得到一个离散的点集。点数趋于无穷多，而欧氏长度趋于零。
32、经无限操作，达到极限时所得到的离散点集称之为Cantor 集。模拟程序function f=cantor(ax,ay,bx,by) c=0.2;d=2; if (bx-ax)c x=ax,bx;y=ay,by;hold on; plot(x,y,LineW idth,5);hold off; cx=ax+( bx-ax)/3; cy=ay-d; dx=bx-(bx-ax)/3; dy=by-d; ay=ay-d; by=by-d; cantor(ax,ay,cx,cy); cantor(dx,dy,bx,by); end Koch （科赫）曲线在一单位长度的线段上对其三等分，将中间段直线换成一
33、个去掉底边的等边三角形，再在每条直线上重复以上操作，如此进行下去直到无穷，就得到分形曲线Koch 曲线。模拟程序function f=Koch(ax,ay,bx,by,c) if (bx-ax)2+(by-ay)2=0&(ex-cx)0)|(alpha=0&(ex-cx)0) alpha=alpha+pi; end dy=cy+sin(alpha+pi/3)*l; dx=cx+cos(alpha+pi/3)*l; Koch(ax,ay,cx,cy,c); Koch(ex,ey,bx,by,c); Koch(cx,cy,dx,dy,c); Koch(dx,dy,ex,ey,c); end 输入：
34、 Koch(0,0,120,0,10); Koch 曲线也可以 “ 随机生成 ” 。在构造的每一步，每次去掉区间中间三分之一的部分，而用与去掉部分构成等边三角形的另两边来代替，再用抛硬币的方法决定新的部分位于被去掉的“ 上边” 或“ 下边” 。经几步以后，得到一个看起来相当不规则的随机Koch 曲线，它仍然保留了Koch 曲线的某些特征，如具有精细的结构，但Koch 曲线具有的严格子相似性已被它所具有的“ 统计自相似性 ” 所取代。模拟程序将 Koch() 中第 13、14 行改为：val=rand(); if val0.95 vv=1; else vv=-1; end dy=cy+sin(a
35、lpha+vv*pi/3)*l; dx=cx+cos(alpha+vv*pi/3)*l; 输入： Koch(0,0,120,0,10); Koch snowflake( 科赫雪花 ) 科赫雪花也称为科赫岛，她的构造是三条科赫曲线的合成，在一个三角形的三条边上，按照科赫曲线的构造过程不断分割替代，就形成科赫雪花。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 13 页，共 14 页 - - - - - - - - - 模拟程序function f=kochsnow(l,c) ax=0;
36、ay=0;bx=l;by=0; cx=l/2; cy=l*sqrt(3)/2; Koch(ax,ay,cx,cy,c); Koch(cx,cy,bx,by,c); Koch(bx,by,ax,ay,c); 输入： kochsnow(120,5); 分形树下面给出树形分形图的MATLAB 模拟程序。该分形树是对称的，利用画线命令line 。模拟程序h=10;d=0.6;th=pi/4;n=9; A=zeros(2,2(n+1)-2);A(:,2)=0;h;C=cos(th) -sin(th);sin(th) cos(th);D=inv(C); for i=2:n B(:,2(i-1)-1:2i-
37、2)=A(:,2(i-1)-1:2i-2); A(:,2i-1:2i+2( i-1)-2)=d*C*B(:,2(i-1)-1:2i-2)+0;h*ones(1,2i-2(i-1); A(:,2i+2(i-1)-1:2(i+1)-2)=d*D*B(:,2(i-1)-1:2i-2)+0;h*ones(1,2i-2(i-1); end for i=1:2n-1 L=line(A(1,2*i-1 2*i),A(2,2*i-1 2*i);set(L,LineWidth,2); end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 14 页，共 14 页 - - - - - - - - -