t=0:0.0001:0.1; %时间间隔为0.0001，说明采样频率为10000Hz

x=square(2*pi*1000*t); %产生基频为1000Hz的方波信号

n=randn(size(t)); %白噪声

f=x+n; %在信号中加入白噪声

figure(1);

subplot(2,1,1);

plot(f); %画出原始信号的波形图

ylabel('幅值(V)');

xlabel('时间(s)');

title('原始信号');

y=fft(f,1000); %对原始信号进行离散傅里叶变换，参加DFT采样点的个数为1000

subplot(2,1,2);

m=abs(y);

f1=(0:length(y)/2-1)'*10000/length(y);%计算变换后不同点对应的幅值

plot(f1,m(1:length(y)/2));

ylabel('幅值的模');

xlabel('时间(s)');

title('原始信号傅里叶变换');

%用周期图法估计功率谱密度

p=y.*conj(y)/1000; %计算功率谱密度

ff=10000*(0:499)/1000; %计算变换后不同点对应的频率值

figure(2);

plot(ff,p(1:500));

ylabel('幅值');

xlabel('频率(Hz)');

title('功率谱密度（周期图法）');

功率谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合matlab，我做一个粗略介绍。功率谱估计可以分为经典谱估计方法与现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计.都可以编程实现，很简单。在matlab中，周期图法可以用函数periodogram实现。

　　但是周期图法估计出的功率谱不够精细，分辨率比较低。因此需要对周期图法进行修正，可以将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数据的功率谱估计。

这2种称为分段平均周期图法，一般后者比前者效果好。加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加一个非矩形窗进行预处理，然后在按分段平均周期图法估计功率谱。相对于分段平均周期图法，加窗平均周期图法可以减小频率泄漏，增加频峰的宽度。welch法就是利用改进的平均周期图法估计估计随机信号的功率谱，它采用信号分段重叠，加窗，FFT等技术来计算功率谱。与周期图法比较，welch法可以改善估计谱曲线的光滑性，大大提高谱估计的分辨率。

　　现代谱估计主要针对经典谱估计分辨率低和方差性不好提出的，可以极大的提高估计的分辨率和平滑性。可以分为参数模型谱估计和非参数模型谱估计。参数模型谱估计有AR模型，MA模型，ARMA模型等；非参数模型谱估计有最小方差法和MUSIC法等。由于涉及的问题太多，这里不再详述，可以参考有关资料。

matlab中，现代谱估计的很多方法都可以实现。music方法用pmusic命令实现；pburg函数利用burg法实现功率谱估计；pyulear函数利用yule-walker算法实现功率谱估计等等。

　　另外，sptool工具箱也具有功率谱估计的功能。窗口化的操作界面很方便，而且有多种方法可以选择。

经典功率谱估计 

直接法：

直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。

Matlab代码示例：

clear;

Fs=1000; %采样频率

n=0:1/Fs:1;

%产生含有噪声的序列

xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));

window=boxcar(length(xn)); %矩形窗

nfft=1024;

[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法

plot(f,10*log10(Pxx));



间接法：

间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。

Matlab代码示例：

clear;

Fs=1000; %采样频率

n=0:1/Fs:1;

%产生含有噪声的序列

xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); 

nfft=1024;

cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数

CXk=fft(cxn,nfft);

Pxx=abs(CXk);

index=0:round(nfft/2-1);

k=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

plot(k,plot_Pxx);



改进的直接法：

对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。

1. Bartlett法

Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。

Matlab代码示例：

clear；

Fs=1000;

n=0:1/Fs:1;

xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));

nfft=1024;

window=boxcar(length(n)); %矩形窗

noverlap=0; %数据无重叠

p=0.9; %置信概率

[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);

index=0:round(nfft/2-1);

k=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

plot_Pxxc=10*log10(Pxxc(index+1));

figure(1)

plot(k,plot_Pxx);

pause;

figure(2)

plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

改进的直接法：

对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。

1. Bartlett法

Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。

Matlab代码示例：

clear；

Fs=1000;

n=0:1/Fs:1;

xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));

nfft=1024;

window=boxcar(length(n)); %矩形窗

noverlap=0; %数据无重叠

p=0.9; %置信概率

[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);

index=0:round(nfft/2-1);

k=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

plot_Pxxc=10*log10(Pxxc(index+1));

figure(1)

plot(k,plot_Pxx);

pause;

figure(2)

plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

2. Welch法

Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。 

Matlab代码示例：

clear;

Fs=1000;

n=0:1/Fs:1;

xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));

nfft=1024;

window=boxcar(100); %矩形窗

window1=hamming(100); %海明窗

window2=blackman(100); %blackman窗

noverlap=20; %数据无重叠

range='half'; %频率间隔为[0 Fs/2]，只计算一半的频率

[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);

[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);

[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);

plot_Pxx=10*log10(Pxx);

plot_Pxx1=10*log10(Pxx1);

plot_Pxx2=10*log10(Pxx2);



figure(1)

plot(f,plot_Pxx);



pause;



figure(2)

plot(f,plot_Pxx1);



pause;



figure(3)

plot(f,plot_Pxx2);

做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。

自相关

x(n)和y(n)为统计的随机序列，-∞其中E{·}为预期的数值操作，且xcorr函数只能估计有限序列。c = xcorr(x,y) 返回矢量长度为2*N-1互相关函数序列，其中x和y的矢量长度均为N，如果x和y的长度不一样，则在短的序列后补零直到两者长度相等

输出矢量C通过c(m) = Rxy(m-N), m=1, ..., 2N-1.得到。

通常，互相关函数需要正规化来得到更准确的估计值。

c = xcorr(x) 为矢量x的自相关估计；

c = xcorr(x,y,'option') 为有正规化选项的互相关计算；其中选项为

"biased"为有偏的互相关函数估计；

"unbiased"为无偏的互相关函数估计；

"coeff"为0延时的正规化序列的自相关计算；

"none"为原始的互相关计算；

若需了解更多的有偏和无偏相关估计，请参阅文献[1].

c= xcorr(x,'option')特指以上某个选项的自相关估计。

c = xcorr(x,y,maxlags) 返回一个延迟范围在[-maxlags,maxlags]的互相关函数序列，输出c的程度为2*maxlags+1.

c = xcorr(x,maxlags) 返回一个延迟范围在[-maxlags,maxlags]的自相关函数序列，输出c的程度为2*maxlags+1.

c = xcorr(x,y,maxlags,'option') 同时指定maxlags和option的互相关计算.

c = xcorr(x,maxlags,'option') 同时指定maxlags和option的自相关计算.

[c,lags] = xcorr(...)返回一个在c进行相关估计的延迟矢量lag，其范围为[-maxlags:maxlags],当maxlags没有指定时，其范围为[-N+1,N-1]

lags就是信号的延时或者超前啊，两个信号的相关性是相对于一定超前和滞后而言的。相关算法就是用移位相乘来体现信号之间的相似度，这里包括幅值和频率。你看看帮助里面的例子：

x = [1,2i,3]; y = [4,5,6];

[c1,lags] = xcorr(x,y)

功率谱密度