离散数字信号的频谱分析程序共享（带注释）
在这个论坛学到了很多东西，通过五一七天的突击，我从一窍不通到能编写小程序了，不过我一直没找到离散数字信号的频谱分析完整程序，为了让新手借鉴，这里贴出我编的程序，请高手指教^_^%fft变换,获得采样数据基本信息，原始曲线图，fft变换频谱图
%这里的向量都是行向量
%注意：采样频率的给定直接会影响频谱图和最终的周期
clear;
close all;
load data.txt %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版
A=data; %将测量数据赋给A，此时A为N×2的数组
x=A(:,1); %将A中的第一列赋值给x，形成时间序列
y=A(:,2); %将A中的第二列赋值给y，形成被测量序列
%注意：这里x、y都是列向量
%显示数据基本信息
fprintf('\n数据基本信息：\n')
fprintf(' 采样点数 = %7.0f \n',length(x)) %输出采样数据个数
fprintf(' 采样时间 = %7.3f s\n',max(x)-min(x)) %输出采样耗时
fprintf(' 采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))%输出采样频率
fprintf(' 最小速度 = %7.3f m/s\n',min(y)) %输出本次采样被测量最小值
fprintf(' 平均速度 = %7.3f m/s\n',mean(y)) %输出本次采样被测量平均值
fprintf(' 速度中值 = %7.3f m/s\n',median(y)) %输出本次采样被测量中值
fprintf(' 最大速度 = %7.3f m/s\n',max(y)) %输出本次采样被测量最大值
fprintf(' 标准方差 = %7.3f \n',std(y)) %输出本次采样数据标准差
fprintf(' 协 方 差 = %7.3f \n',cov(y)) %输出本次采样数据协方差
fprintf(' 自相关系数 = %7.3f \n\n',corrcoef(y)) %输出本次采样数据自相关系数
%显示原始数据曲线图
figure;
subplot(2,1,1)
plot(x,y) %显示原始数据曲线图
axis() %设置坐标最大值，优化作用
xlabel('时间 (s)')
ylabel('被测量 (单位)')
title('原始信号')
%傅立叶变换
y=y-mean(y); %消去直流分量，使频谱更能体现有效信息
Fs=2000; %得到原始数据data.txt时，仪器的采样频率
N=10000; %data.txt中的被测量个数，即采样个数
z=fft(y); %注意因为y是列向量，所以这里z是列向量
%功率谱分析
f=(0:N-1)'*Fs/N; %因为z是列向量，所以这里f要转换成列向量
Mag=2*abs(z)/N; %幅值
Pyy=Mag.^2; %功率谱
%显示功率谱－频率图
subplot(2,1,2)
plot(f(1:N/2),Pyy(1:N/2),'r') %显示功率谱－频率
% |
% 将这里的Pyy改成Mag就是 幅值-频率图 了
axis()%坐标优化，可有可无
xlabel('频率 (Hz)')
ylabel('功率谱')
title('功率谱 - 频率')
grid on
%返回最大功率谱对应的频率和周期值
=max(Pyy(1:N/2));
fprintf('\n傅立叶变换结果：\n')
fprintf(' FFT_f = %1.3f Hz\n',f(b)) %输出最大值对应的频率
fprintf(' FFT_T = %1.3f s\n',1/f(b)) %输出最大值对应的周期
[ 本帖最后由 wufashengcun 于 2007-5-7 23:39 编辑 ]
上面程序的改进版
呵呵，发出来这么久，怎么都没人对我的帖子提出疑义呢？边学习边自己改进了，下面是改进后的程序：
%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数据基本信息：\n')
fprintf(' 采样点数 = %7.0f \n',length(x)) %输出采样数据个数
fprintf(' 采样时间 = %7.3f s\n',max(x)-min(x)) %输出采样耗时
fprintf(' 采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))%输出采样频率
fprintf(' 最小速度 = %7.3f m/s\n',min(y)) %输出本次采样被测量最小值
fprintf(' 平均速度 = %7.3f m/s\n',mean(y)) %输出本次采样被测量平均值
fprintf(' 速度中值 = %7.3f m/s\n',median(y)) %输出本次采样被测量中值
fprintf(' 最大速度 = %7.3f m/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() %优化坐标，可有可无
xlabel('时间 (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);
%频谱分析
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()
xlabel('频率 (Hz)')
ylabel('能量')
title('频谱图(频域)')
grid on;
%返回最大能量对应的频率和周期值
=max(Pyy(1:N/2));
fprintf('\n傅立叶变换结果：\n')
fprintf(' FFT_f = %1.3f Hz\n',f(b)) %输出最大值对应的频率
fprintf(' FFT_T = %1.3f s\n',1/f(b)) %输出最大值对应的周期
[ 本帖最后由 wufashengcun 于 2007-6-4 13:03 编辑 ] 不错，正是我想要找的东西！ 我也在学，有帮助啊:@)
很有启发，但有连续信号的吗？
同为初学者，没你进步快，惭愧。
页:
[1]