目录
【若觉文章质量良好且有用，请别忘了点赞收藏加关注，这将是我继续分享的动力，万分感谢！】
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. 场景仿真：
输入数据：
得到的能谱
得到的三分之一倍频程能谱
【若觉文章质量良好且有用，请别忘了点赞收藏加关注，这将是我继续分享的动力，万分感谢！】