注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

神瑛侍者的博客

 
 
 

日志

 
 

【转载】我的Matlab幅值谱和功率谱计算函数  

2013-05-08 21:31:07|  分类: 默认分类 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

function [A, P, f] = myfft(y, fs, unit)
% y 必须为向量
% fs 为采样频率
% unit 为数据y的单位
% f 频谱图的频率变量
% A 为对应于 f 处的幅值
% P 为对应于 f 处的功率
% 例如对采样频率为1000Hz的振动数据y(单位为mm)进行傅里叶变换,调用格式为
% [A, P, f] = myfft(y, 1000, 'mm')
%
len = length(y);
for i = 1:20
    if (len < 2^i)      % 如果2^i大于y的长度
        N = 2^(i-1);    % 则截取y的前2^(i-1)个数据参与计算,后面的部分忽略掉
        break;
    end
end

t = linspace(0, len/fs, len);   % 画时域波形
figure
plot(t, y)
title('时域波形')
xlabel('时间 t /s')
ylabel(['幅值 A /',unit])
grid on

Y = fft(y, N);          % 快速傅里叶变换Y = Y(k),k = 1,2,...N
f = fs/N*(0:N/2);       % 频率步长Δf = 1/T = 1/((1/fs)*N) = fs/N

Ay = abs(Y)/N;     % 双边幅值谱
A = [Ay(1), 2*Ay(2:N/2+1)]; % 单边幅值谱
%A = 2*Ay(1:N/2+1);
figure
plot(f, A)
title('幅值谱')
xlabel('频率 f /Hz')
ylabel(['幅值 A /',unit])
grid on

Py = abs(Y).^2/N^2;     % 双边功率谱
P = [Py(1), 2*Py(2:N/2+1)];      % 单边功率谱
%P = 2*Py(1:N/2+1);
figure
plot(f, P)
title('功率谱')
xlabel('频率 f /Hz')
ylabel(['功率 P /',unit,'^2'])
grid on

  评论这张
 
阅读(38)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017