离散信号MATLAB频谱分析程序

离散信号MATLAB频谱分析程序【ZZ】%FFT变换,获得采样数据基本信息，时域图，频域图%这里的向量都用行向量，假设被测变量是速度，单位为m/sclear;closeall;loaddata.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数据基本信息：\n')fprintf('采样点数=%7.0f\n',length(x))%输出采样数据个数fprintf('采样时间=%7.3fs\n',max(x)-min(x))%输出采样耗时fprintf('采样频率=%7.1fHz\n',length(x)/(max(x)-min(x)))%输出采样频率fprintf('最小速度=%7.3fm/s\n',min(y))%输出本次采样被测量最小值fprintf('平均速度=%7.3fm/s\n',mean(y))%输出本次采样被测量平均值fprintf('速度中值=%7.3fm/s\n',median(y))%输出本次采样被测量中值fprintf('最大速度=%7.3fm/s\n',max(y))%输出本次采样被测量最大值fprintf('标准方差=%7.3f\n',std(y))%输出本次采样数据标准差fprintf('协方差=%7.3f\n',cov(y))%输出本次采样数据协方差fprintf('自相关系数=%7.3f\n\n',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('时间(s)');ylabel('被测变量y');title('原始信号(时域)');gridon;%傅立叶变换y=y-mean(y);%消去直流分量，使频谱更能体现有效信息Fs=2000;%得到原始数据data.txt时，仪器的采样频率。其实就是length(x)/(max(x)-min(x));N=10000;%data.txt中的被测量个数，即采样个数。其实就是length(y);z=fft(y);%频谱分析f=(0:N-1)*Fs/N;Mag=2*abs(z)/N;%幅值，单位同被测变量yPyy=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)))])xlabel('频率(Hz)')ylabel('能量')title('频谱图(频域)')gridon;%返回最大能量对应的频率和周期值[ab]=max(Pyy(1:N/2));fprintf('\n傅立叶变换结果：\n')fprintf('FFT_f=%1.3fHz\n',f(b))%输出最大值对应的频率fprintf('FFT_T=%1.3fs\n',1/f(b))%输出最大值对应的周期*************************************************************************%%FFT实践及频谱分析%%*************************************************************************%%*************************************************************************%%***************1.正弦波****************%fs=100;%设定采样频率N=128;n=0:N-1;t=n/fs;f0=10;%设定正弦信号频率%生成正弦信号x=sin(2*pi*f0*t);figure(1);subplot(231);plot(t,x);%作正弦信号的时域波形xlabel('t');ylabel('y');title('正弦信号y=2*pi*10t时域波形');grid;%进行FFT变换并做频谱图y=fft(x,N);%进行fft变换mag=abs(y);%求幅值f=(0:length(y)-1)'*fs/length(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);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;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');ylabel('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('均方根谱');title('矩形波均方根谱');grid;%求功率谱power=sq.^2;figure(2);subplot(234);plot(f,power);xlabel('频率(Hz)');ylabel('功率谱');title('矩形波功率谱');grid;%求对数谱ln=log(sq);figure(2);subplot(235);plot(f,ln);xlabel('频率(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);plot(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)');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(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一个计算图像分形维数的软件这个是计算图象的分形维数的软件,可在WINDOWS下安装运行,运行时,新建一个project,在其中导入图象即可下载地址：=596还有萝卜驿站中也有很多不错的东西：里面有一个盒维的MATLAB程序一个分形的网站：分形频道2003看下里面的留言区，有些东西哦一个求盒维数的Matlab程序：%本程序只能计算时间序列的盒维数%box11cleardata1;%存放时间序列的数据文件，共6列c=0;n=2^10;ll=1:1383-1220+1;b(:,1)=ll';b(:,2)=aaa(1220:1383,6)*100;%aaa是data1中的函数，里面是6列时间序列%划分不同大小的格子，并统计不同量级格子中含有数据的格子数k=2;forp=1:9...k=2^p...stepx=(max(b(:,1))-min(b(:,1))+1)/k;...stepy=(max(b(:,2))-min(b(:,2))+1)/k;...mm(p)=0;...fori=1:k......i......forj=1:k.........lowx=stepx*(i-1)+min(