参考文章：百度文库文章

但是该文章有很多部分，不知是matlab更新了还是咋地，有错误，踩了很多坑，仅做参考。
产生高斯白噪声的方法在“信号处理”专栏中有写，其时域和频域图如下：



用matlab中求自（互）相关的xcorr函数，参考帮助文档
[r,lags]=xcorr(z,'biased');

其中z是上述高斯白噪声，r是自相关函数，lags是时间偏移量（索引），尤其注意’biased’参数，这是调试了半天才发现的问题。

帮助文档里描述如下：



相当于算出来的r除以了样本数，这是离散自相关函数计算步骤。

如果不加这个参数，自相关函数值会变得很大。

算出自相关函数之后，根据维纳辛钦定理，信号的功率谱密度和信号的自相关函数是一个傅里叶变换对。使用fft函数求自相关函数的傅里叶变换：
pdv=fft(r);
pdv=abs(fftshift(pdv))./length(pdv);

fftshift是频谱矫正函数，将fft后的结果以中心为分界，两边分别做镜像，得出来的结果再除以样本数，即双边频谱图。

而这个频谱图就是原信号的功率谱密度。
plot((0:length(pdv)-1)*L/length(pdv)-L/2,pdv)

注意这里横坐标是频率，需要做映射处理，L为样本数量


高斯白噪声其实看不出来啥，下面用一个没加噪声的纯信号来展示：
y=12*cos((2*pi)*100.*t)+15*cos((2*pi)*150.*t)+18*cos((2*pi)*210.*t)+10;




上面的功率谱是由双边频谱的平方（除直流）计算出来的，下面的是用维纳辛钦定理计算出来的，可以看到功率谱几乎相同。

信号的自相关函数如下：



编者水平有限，很多知识点都是似懂非懂，如有错误欢迎指出！

附代码：
fs=1000;%采样频率hz
T_N=1.5;%总时间s
t=1/fs:1/fs:T_N;%时间向量
L=T_N*fs;%样本长度
y=12*cos((2*pi)*100.*t)+15*cos((2*pi)*150.*t)+18*cos((2*pi)*210.*t)+10;%信号
subplot(2,1,1);
plot(t,y);
xlabel("时间/s")
ylabel("幅度/v")
title("时域")
fft_y=fft(y);%快速傅里叶变换
P = abs(fft_y/L);%取幅频特性，除以L
P = P(1:L/2+1);%截取前半段
P(2:end-1)=2*P(2:end-1);%单侧频谱非直流分量记得乘以2
f = fs*(0:(L/2))/L;%频率，最多到一半（奈奎斯特采样定理）
subplot(2,1,2);
plot(f,P);
xlabel("频率/Hz")
ylabel("幅度/v")
title("单边频谱")

figure(4)
[r,lags]=xcorr(y,'biased');%得到自相关函数的幅度和偏移量
subplot(2,1,1)
plot(lags,r)
xlabel("时间偏移/s")
ylabel("相关程度")
title("（自相关函数）")

fft_y=fftshift(fft_y);%频谱矫正
powerpu=abs(fft_y./L).^2;%双边功率谱
% subplot(2,1,1)
% plot((1:size(powerpu,2))-751,powerpu)
% xlabel("频率/Hz")
% ylabel("功率/W")
% title("功率谱")

pdv=fft(r,size(r,2));%对自相关函数快速傅里叶变换
pdv=abs(fftshift(pdv));%频谱矫正，让正半轴部分和负半轴部分的图像分别关于各自的中心对称，得到双边谱
subplot(2,1,2)
plot((1:length(pdv))*L/length(pdv)-L/2,pdv./length(powerpu))
xlabel("频率/Hz")
ylabel("功率谱W/Hz")
title("（功率谱）")

% z1=0.1*randn(1,201);%产生方差N(0,0.12)高斯白噪声
% [r1,lags]=xcorr(z1,'unbiased');%自相关函数的估计
% plot(lags,r1);
% f1=fft(r1);
% f2=fftshift(f1);%频谱校正
% l1=((0:length(f2)-1)*200/length(f2)-100);%功率谱密度x轴
% y4=abs(f2);
% figure(2)
% plot(l1,y4);

总述

一、区分信号类型
根据式$E=\int_{-\infty }^{\infty }{{{s}^{2}}\left( t \right)}dt$计算信号能量（作用在单位电阻上的电压信号 释放的能量）可以将信号分为：

功率信号：能量无限，不能用能量表示，所以用平均功率表示；
能量信号：能量有限，平均功率为0；

二、功率信号的分析

频谱（离散）：$C\left( n{{f}_{0}} \right)=\frac{1}{{{T}_{0}}}\int_{-T/2}^{T/2}{s\left( t \right){{e}^{-j2\pi n{{f}_{0}}t}}dt}$（单位：V）

含义：

周期功率信号幅值（频率为${{f}_{0}}$ ）经过傅里叶级数展开，被多个离散倍频$n{{f}_{0}}$表征，各频点的幅值$C\left( n{{f}_{0}} \right)$也即该频点的贡献权系数。

参考网上一张经典的解释图：


功率谱密度（连续）：$P\left( f \right)={{\sum{\left| C\left( f \right) \right|}}^{2}}\delta \left( f-n{{f}_{0}} \right),C\left( f \right)=\left\{ \begin{matrix}
C\left( n{{f}_{0}} \right),f=n{{f}_{0}}  \\
0,  \\
\end{matrix} \right.$

含义：


将信号的功率按照频点贡献铺在频谱之上；
因其能量是无穷的，所以不能把能量铺上去，只能用有限的功率；
对功率谱密度进行积分，能得到局部频段承载的功率；
相比功率信号的频谱突出各频点对功率信号的信号幅值的贡献，功率谱密度突出各频点对功率信号的功率的贡献。

三、能量信号的分析

频谱密度（连续）：$S\left( f \right)=\int_{-\infty }^{\infty }{s\left( t \right){{e}^{-j2\pi ft}}dt}$（单位：V/Hz）

含义：


通过傅里叶变换将能量信号转换到连续频域上；
但因能量有限，不能使用离散贡献频点权系数（几乎为0），只能使用频谱密度来表征。


能量谱密度：$G\left( f \right)={{\left| S\left( f \right) \right|}^{2}}$（单位：J/Hz）

含义：


将信号能量铺在频谱之上；
对能量谱密度进行局部积分，能得到局部频段承载的能量；
相比能量信号的频谱密度突出连续频点对功率信号的信号幅值的贡献，能量谱密度突出连续频点对能量信号的能量的贡献。

四、相关函数

自相关函数  ：$R\left( \tau  \right)=\int_{-\infty }^{\infty }{s\left( t \right)s\left( t+\tau  \right)}dt$

含义：表征信号自己与自己的相似关系，根据公式也能看出，原信号不变，拿着一个复制信号先进行时移再乘积积分（卷积操作），就是在求信号内时间维的相似性，$\tau =0$ 时肯定最相似，周期信号肯定有很多相似点（自相关函数峰值）。

特点：


互相关函数：${{R}_{1,2}}\left( \tau  \right)=\int_{-\infty }^{\infty }{{{s}_{1}}\left( t \right){{s}_{2}}\left( t+\tau  \right)}dt$

含义：表征两信号的相似关系

特点也同上表

注：上文具体公式推导见樊昌信的《通信原理》第六版P17-33

	

来源：声振测试微信公众号，作者：于长帅。
频率响应函数表征了测试系统对给定频率下的稳态输出与输入的关系，是试验分析和仿真分析的基础。频响函数概念及分析流程
互谱的一个重要应用是计算线性系统的频率响应函数。设x(t ) 为某点输入的平稳随机信号，y(t ) 为任一点的响应，也是随机平稳信号，则振动系统的频率响应函数的互功率谱密度函数除以自功率谱密度函数的商，即：
式中，Sxx(f ) 和Sxy(f ) 分别为随机振动信号的自功率谱密度函数和激励与响应信号的互功率谱密度函数的估计，无论系统是稳态的、非稳态的，还是确定性的，上式计算频响函数都是唯一正确的。
频率响应函数是复函数，它是被测系统的动力特性在频域内的表现形式，也就是被测系统本身对输入信号在频域中的传递特性的描述。输入信号的各频率成分通过该系统时，频响函数对它们的一些频率成分进行了放大，对另一些频率成分进行了衰减，经过加工后得到输出信号的新频率成分的分布，因此频响函数对结构的动力特性测试具有特殊重要的意义。
图1 频响函数分析流程
自相关函数：
为了反映信号自身取值随自变量时间前后变化的相似性对于信号x(t )而言，其自相关函数可以表示为：
式中，τ 为自相关函数的时延量；T 为信号周期。对信号进行自相关分析，就是对信号x(t )延迟时间τ 后的信号x(t+τ)，然后对x(t )和x(t+τ)做卷积计算，所得结果即为x(t )的自相关函数。
互相关函数：
际中将信号y(t ) 移动时间τ 得到y(t+τ)，然后再计算x(t )和y(t+τ)的相关性。互相关函数可以写成可写成：
式中，T 为信号x(t ) 和y(t ) 的观测时间，τ 为信号的滞后时间，R(τ)是τ 的函数。
振动信号自相关功率谱密度函数：
自功率谱描述了信号的频率结构，反映了振动能量在各个频率上的分布情况。自功率谱密度函数的定义是自相关的傅立叶变换，如下式：
式中，为振动信号x(t ) 的自相关函数，是时域中的统计量。
振动信号互相关功率谱密度函数：
同理，互功率谱密度函数的定义是互相关函数的傅立叶变换，如下式：
频响函数分析实例
下图用激振器对某产品进行模态测试，产品上布置着加速度传感器，激振器力传感器的时域信号和产品某测点的加速度信号如下图所示。
图2 激振器模态测试
(a) 激振器力时域信号；(b) 某测点加速度信号
图3 激振器力传感器的时域信号和某测点的加速度信号
对该测点的加速度信号进行自相关分析，和对该测点的加速度信号和力信号的互相关分析，分析结果如下图所示。
(a) 加速度测点自相关函数曲线；(b) 力与加速度互相关函数曲线
图4 自相关函数和互相关函数曲线
对自相关函数和互相关函数分别进行傅里叶变化，分析结果如下图所示。
(a)自相关函数傅里叶分析曲线；(b) 互相关函数傅里叶分析曲线
图5 自相关函数和互相关函数进行傅里叶变化
经过以上分析，激振器力信号和加速度测点信号的频响函数曲线如下图所示。
图6 频率响应函数曲线
附录MATLAB程序：
clear
clc
x=xlsread('C:\Users\Administrator\Desktop\FRF_matlab\ACC.xlsx');
%导入要分析的加速度信号
x1=x(1:1500,2);
[b,a]=xcorr(x1,'unbiased');
%对加速度信号进行自相关分析
y=xlsread('C:\Users\Administrator\Desktop\FRF_matlab\force.xlsx');
%导入要分析的输入力信号
y1=y(1:1500,2);
[c,d]=xcorr(x1,y1,'unbiased');
%做互相关分析
plot(c,d);
xlabel('时间/s');
ylabel('加速度/g');
L=1500;
Fs=512;
NFFT=2^nextpow2(L);
Y=fft(b,NFFT)/L;
f=(Fs/2)*linspace(0,1,NFFT/2+1);
Y1=fft(c,NFFT)/L;
%对自相关函数进行傅里叶变化
f=(Fs/2)*linspace(0,1,NFFT/2+1);
%对互相关函数进行傅里叶变化
Z1=2*abs(Y1(1:NFFT/2+1));
Z2=2*abs(Y(1:NFFT/2+1));
loglog(f,Z1./Z2)
%绘制频响函数曲线
xlabel('频率/Hz');
ylabel('加速度/g');
关联阅读：支撑方式对于频率响应函数测量的影响
传递函数、频率响应函数及参数识别详解
锤击法测试中频率响应函数的问题探讨模态频率响应函数的互易性如何理解它
声明：本微信转载文章出于非商业性的教育和科研目的，并不意味着支持其观点或证实其内容的真实性。版权归原作者所有，如转载稿涉及版权等问题，请立即联系我们，我们会予以更改或删除相关文章，保证您的权利！

	

来源：声振测试微信公众号，作者：于长帅。
频率响应函数表征了测试系统对给定频率下的稳态输出与输入的关系，是试验分析和仿真分析的基础。频响函数概念及分析流程
互谱的一个重要应用是计算线性系统的频率响应函数。设x(t ) 为某点输入的平稳随机信号，y(t ) 为任一点的响应，也是随机平稳信号，则振动系统的频率响应函数的互功率谱密度函数除以自功率谱密度函数的商，即：
式中，Sxx(f ) 和Sxy(f ) 分别为随机振动信号的自功率谱密度函数和激励与响应信号的互功率谱密度函数的估计，无论系统是稳态的、非稳态的，还是确定性的，上式计算频响函数都是唯一正确的。
频率响应函数是复函数，它是被测系统的动力特性在频域内的表现形式，也就是被测系统本身对输入信号在频域中的传递特性的描述。输入信号的各频率成分通过该系统时，频响函数对它们的一些频率成分进行了放大，对另一些频率成分进行了衰减，经过加工后得到输出信号的新频率成分的分布，因此频响函数对结构的动力特性测试具有特殊重要的意义。
图1 频响函数分析流程
自相关函数：
为了反映信号自身取值随自变量时间前后变化的相似性对于信号x(t )而言，其自相关函数可以表示为：
式中，τ 为自相关函数的时延量；T 为信号周期。对信号进行自相关分析，就是对信号x(t )延迟时间τ 后的信号x(t+τ)，然后对x(t )和x(t+τ)做卷积计算，所得结果即为x(t )的自相关函数。
互相关函数：
际中将信号y(t ) 移动时间τ 得到y(t+τ)，然后再计算x(t )和y(t+τ)的相关性。互相关函数可以写成可写成：
式中，T 为信号x(t ) 和y(t ) 的观测时间，τ 为信号的滞后时间，R(τ)是τ 的函数。
振动信号自相关功率谱密度函数：
自功率谱描述了信号的频率结构，反映了振动能量在各个频率上的分布情况。自功率谱密度函数的定义是自相关的傅立叶变换，如下式：
式中，为振动信号x(t ) 的自相关函数，是时域中的统计量。
振动信号互相关功率谱密度函数：
同理，互功率谱密度函数的定义是互相关函数的傅立叶变换，如下式：
频响函数分析实例
下图用激振器对某产品进行模态测试，产品上布置着加速度传感器，激振器力传感器的时域信号和产品某测点的加速度信号如下图所示。
图2 激振器模态测试
(a) 激振器力时域信号；(b) 某测点加速度信号
图3 激振器力传感器的时域信号和某测点的加速度信号
对该测点的加速度信号进行自相关分析，和对该测点的加速度信号和力信号的互相关分析，分析结果如下图所示。
(a) 加速度测点自相关函数曲线；(b) 力与加速度互相关函数曲线
图4 自相关函数和互相关函数曲线
对自相关函数和互相关函数分别进行傅里叶变化，分析结果如下图所示。
(a)自相关函数傅里叶分析曲线；(b) 互相关函数傅里叶分析曲线
图5 自相关函数和互相关函数进行傅里叶变化
经过以上分析，激振器力信号和加速度测点信号的频响函数曲线如下图所示。
图6 频率响应函数曲线
附录MATLAB程序：
clear
clc
x=xlsread('C:\Users\Administrator\Desktop\FRF_matlab\ACC.xlsx');
%导入要分析的加速度信号
x1=x(1:1500,2);
[b,a]=xcorr(x1,'unbiased');
%对加速度信号进行自相关分析
y=xlsread('C:\Users\Administrator\Desktop\FRF_matlab\force.xlsx');
%导入要分析的输入力信号
y1=y(1:1500,2);
[c,d]=xcorr(x1,y1,'unbiased');
%做互相关分析
plot(c,d);
xlabel('时间/s');
ylabel('加速度/g');
L=1500;
Fs=512;
NFFT=2^nextpow2(L);
Y=fft(b,NFFT)/L;
f=(Fs/2)*linspace(0,1,NFFT/2+1);
Y1=fft(c,NFFT)/L;
%对自相关函数进行傅里叶变化
f=(Fs/2)*linspace(0,1,NFFT/2+1);
%对互相关函数进行傅里叶变化
Z1=2*abs(Y1(1:NFFT/2+1));
Z2=2*abs(Y(1:NFFT/2+1));
loglog(f,Z1./Z2)
%绘制频响函数曲线
xlabel('频率/Hz');
ylabel('加速度/g');
关联阅读：支撑方式对于频率响应函数测量的影响
传递函数、频率响应函数及参数识别详解
锤击法测试中频率响应函数的问题探讨模态频率响应函数的互易性如何理解它
声明：本微信转载文章出于非商业性的教育和科研目的，并不意味着支持其观点或证实其内容的真实性。版权归原作者所有，如转载稿涉及版权等问题，请立即联系我们，我们会予以更改或删除相关文章，保证您的权利！