http://blog.sciencenet.cn/blog-2777453-1000182.html

%--------------------------------------------------------------------------

%功能：                  谱分析代码

%                    何成文         2016年9月1日

%--------------------------------------------------------------------------

clc;clear;close all;

load('WPM10000.txt');

data=WPM10000;

N=length(data);

figure;plot([1:N],data,'r');title('原始数据');grid on;xlabel('T/s');

%--------------------------------------------------------------------------

F=fft(data,N);

R=abs(F)/N;

Angle=angle(F);

n=[0:1:N-1];

Hz=n/(N*1);

figure

subplot(1,2,1);plot(Hz(1:N/2),R(1:N/2),'r');title('傅里叶幅度谱');xlabel('X/Hz');ylabel('幅度');grid on;

subplot(1,2,2);plot(Hz(1:N/2),Angle(1:N/2),'b');title('傅里叶相位谱');xlabel('X/Hz');ylabel('相位');grid on;

%--------------------------------------------------------------------------

Y=fft(data,N);                                                             %傅里叶变换！

Y1=(abs(Y)).^2/N;                                                      %Y1为周期图谱计算公式！

n=[0:1:N-1];

Hz=n/(N*1);                                                                %1是抽样时间间隔，本文默认为1！

figure;plot(Hz(1:N/2),Y1(1:N/2),'b');                         %表现形式1

xlabel('X/Hz');ylabel('幅度');title('WPM噪声信号的傅里叶周期谱图');grid on;

figure;plot(log(Hz(1:N/2)),log(Y1(1:N/2)),'r');           %表现形式2-取对数！

xlabel('log(f/Fs)');ylabel('log(PSD)');title('WPM噪声信号的傅里叶周期谱图');grid on;

%--------------------------------------------------------------------------

来源：http://blog.sciencenet.cn/blog-2777453-1000182.html