声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1569|回复: 3

[FFT] 贴出我的频谱程序和结果,出来的结果不圆滑,有很多折角。

[复制链接]
发表于 2010-1-12 14:11 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
function myfrequencyplot
% all the files are save as t and s
clc;
clear;

load result_6000_30.mat;
length(t)
%%%-----------------------Program 1----------------------------------
%FFT变换,获得采样数据基本信息,时域图,频域图
%这里的向量都用行向量,假设被测变量是速度,单位为m/s
x = linspace(t(1),t(end),1e2);% 时间序列
x=x';                                            %将列向量变成行向量
y = interp1(t,s(:,4),x);
                                    %将s中的第4列赋值给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))                   %输出本次采样数据自相关系数
  
%傅立叶变换
y=y-mean(y);                      %消去直流分量,使频谱更能体现有效信息
Fs=length(x)/(max(x)-min(x));     %得到原始数据data.txt时,仪器的采样频率。其实就是length(x)/(max(x)-min(x));
N=length(y);                      %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,故这里有很多表达方式

%显示频谱图(频域)

plot(f(1:N/2),Mag(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('频谱图(频域)')
grid on;

%返回最大能量对应的频率和周期值
[a b]=max(Mag(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))             %输出最大值对应的周期
频谱图.jpg
回复
分享到:

使用道具 举报

发表于 2010-1-12 14:25 | 显示全部楼层

回复 楼主 胡晓宇 的帖子

你处理的是什么信号,周期的还是随机的,参见:
http://forum.vibunion.com/forum/thread-89702-1-1.html
 楼主| 发表于 2010-1-12 14:56 | 显示全部楼层

回复 沙发 hcharlie 的帖子

非常您的及时答复。由于以前没有接触过信号处理方面的内容,这部分数据处理已将难为我快两周时间了。
我处理的是转子动力学的结果。应该属于随机信号吧,因为转子的响应并不是周期的。
但是我保存结果仍然是按照转子转动的周期来保存的,结果总共保存了50个转子转动周期时间的数据,不知道这样的样本量够不够,由于采用的是变步长积分,故对结果又进行了插值。
发表于 2010-1-14 22:38 | 显示全部楼层
y=y-mean(y);                      %消去直流分量,使频谱更能体现有效信息
上面这条语句,对频谱是没有影响的吧?因为它相当于信号的所有点减去一个常数,而信号变化快慢是没有改变的,所以对频率也就是得到的频谱图是没有影响的,不知道我说的对不对?
希望讨论。

Mag=2*abs(z)/N;                 %幅值,单位同被测变量y
这句是什么意思?如果不除以N,你再观察一下频谱图。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-14 14:15 , Processed in 0.072966 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表