基于MATLAB的频谱、能量谱、三分之一倍频程分析

目录

1. 相关概念介绍

2. 代码实现

3. 场景仿真：

---

【若觉文章质量良好且有用，请别忘了点赞收藏加关注，这将是我继续分享的动力，万分感谢！】

1. 相关概念介绍

以信号为例，信号在时域下的图形可以显示信号如何随着时间变化，而信号在频域下的图形（一般称为频谱）可以显示信号分布在哪些频率及其比例。频域的表示法除了有各个频率下的大小外，也会有各个频率的相位，利用大小及相位的资讯可以将各频率的弦波给予不同的大小及相位，相加以后可以还原成原始的信号。

时域：描述数学函数或物理信号对时间的关系；例如一个信号的时域波形可以表达信号随着时间的变化[1]。这符合现实世界中人们对信号的认识。

频域：是指在对函数或信号进行分析时，分析其和频率有关部分，而不是和时间有关的部分，和时域一词相对[1]。这不符合现实世界中人们对信号的感官认识，但其为波的形式更符合现实世界中信号的存在，可以辅助人们完成对现实世界中信号的处理，由此产生更多的概率，如频谱、能谱、功率谱、倍频程谱等。

频谱：一个信号是由哪些频率的弦波所组成，也可以看出各频率弦波的大小及相位等信息。一般通过傅里叶变换将时域信号处理成频域信号，获得各个正弦信号的幅度和相位。所得的结果会是分别以幅度及相位为纵轴，频率为横轴的两张图，不过有时也会省略相位的信息，只有不同频率下对应幅度的资料。有时也以“幅度频谱”表示幅度随频率变化的情形，“相位频谱”表示相位随频率变化的情形。

其频谱、能谱、功率谱、倍频程谱的相关概念可看这篇博客：频谱、能谱、功率谱、倍频程谱、1/3 倍频程谱_liyuanbhu的博客-CSDN博客_1/3倍频程

2. 代码实现

这里提供能谱和三分之一倍频程能谱相应的代码（代码不易，请点赞+收藏）：

clc
clear
close all
%% 实验
% xn(:,1) = 1:320000;
% fs = 32000;      % 采样频率
% ts = 1/fs;     % 取样时间
% L = size(xn,1);     % 总取样次数
% t = (0:L-1) * ts;
% xn(:,2) = sin(2*pi*1000*t) + sin(2*pi*5000*t);

%% 正文
xn = xlsread('模拟数据.csv');

% 第一列为采样时刻，第二列为采样数值

%%%%%%%%%%%%%%%%%% 快速傅里叶变换 %%%%%%%%%%%%%%%%%%%%%

fs = 32000;      % 采样频率
ts = 1/fs;     % 取样时间
L = size(xn,1);     % 总取样次数
t = (0:L-1) * ts;
x = xn(:,2);
x = detrend(x,0);

%% 傅里叶变换
nfft = 2^15;

Y1 = fft(x,nfft);     % 傅里叶变换

f1 = (0:nfft-1)/(nfft-1)*fs;      % 频率向量

f = f1(1:nfft/2);       % 根据对称性取半

Y = Y1(1:nfft/2)*2/nfft;   % 根据对称性取半

YE = abs(Y); % 频域中的能量

ji_Ya = 2*10^(-5);  % 声压级基底

%% 计算频谱声压级
Ya_1 = 20 * log10(YE/ji_Ya);  % 计算频谱声压级
    
%% 计算总声压级
% 找出10-8k
YE_z = YE(find(f>=10&f<8000));
Z_Ya = 20 * log10(sqrt(sum(YE_z.^2))/ji_Ya)-3

%% 三分之一倍频程
% 定义三分之一倍频程的中心频率fc
fc = [20 25 31.5 40 50 63 80 100 125 160 200 ...
    250 315 400 500 630 800 1000 1250 1600 2000 ...
    2500 3150 4000 5000 6300 8000 10000 12500 16000];
% 下限频率
fl = round(fc/(2^(1/6)));
% 上限频率
fu = round(fc*(2^(1/6)));

fu(end) = f(end);   % 修复fu，末尾变为16000

%% 
% 频率向量f中有L/2个数据，对应的频率是(0：L/2-1)/(L/2-1)*fs/2;
nl = round(fl*2/fs*(nfft/2-1) + 1); % 下限频率对应的频率向量的序号
nu = round(fu*2/fs*(nfft/2-1) + 1); % 上限频率对应的频率向量的序号
nc = length(fc); % 中心频率的长度


for i = 1:nc
    nn = zeros(1,nfft);
    nn(nl(i):nu(i)) = Y1(nl(i):nu(i));
    nn(end-nu(i)+1:end-nl(i)+1) = Y1(end-nu(i)+1:end-nl(i)+1);
    cc = ifft(nn);
    YE_C(i) = sqrt(var(real(cc(1:nfft))));     % 求取1/3倍频程
%     YE_C(i) = sqrt(sum(YE(nl(i):nu(i)).^2)/2);    % 求取第i个中心频率的能量：频带的平均能量；
end

Ya_2 = 20 * log10(YE_C/ji_Ya); % 计算中心频率的声压级


%% 画图
a(1) = figure(1);
set(gca,'FontSize',10);
plot(0:length(t)-1,x);
title('原始信号');
xlabel('采样序号/采样频率为32000Hz');
ylabel('x');
set(gcf,'position',[100,100, 700, 400]); %设定figure的位置和大小 get current figure
set(gcf,'color','white'); %设定figure的背景颜色

a(2) = figure(2);
plot(f,Ya_1); 
set(gca,'FontSize',10);
title('频谱');
xlabel('频率/Hz');
ylabel('声压级/Db');
set(gcf,'position',[100,100, 700, 400]); %设定figure的位置和大小 get current figure
set(gcf,'color','white'); %设定figure的背景颜色

a(3) = figure(3);
set(gcf,'color','white'); %设定figure的背景颜色
set(gca,'FontSize',10);
plot(fc,Ya_2); 
title('1/3倍频程');
xlabel('中心频率/Hz');
ylabel('声压级/Db');
set(gcf,'position',[100,100, 700, 400]); %设定figure的位置和大小 get current figure



%% 导出数据
write_all = zeros(L,6);
write_all(:,1:2) = xn;
write_all(1:nfft/2,3) = f;
write_all(1:nfft/2,4) = Ya_1;
write_all(1:nc,5) = fc;
write_all(1:nc,6) = Ya_2;

xlswrite('matlab实验结果.xlsx',write_all);

saveas(a(1),'原始信号.bmp');
saveas(a(2),'频谱.bmp');
saveas(a(3),'三分之一倍频程.bmp');

3. 场景仿真：

输入数据：

得到的能谱 

 得到的三分之一倍频程能谱

【若觉文章质量良好且有用，请别忘了点赞收藏加关注，这将是我继续分享的动力，万分感谢！】