三分之一倍频程谱

三分之一倍频程谱是一种频率分析方法，它具有谱线少频带宽的特点。

倍频程实际上是频域分析中频率的一种相对尺度。倍频程谱是由一系列频率点以及对应这些频率点附近的频带内信号的平均幅值（有效值）所构成。这些频率点称为中心频率fc，中心频率附近的频带出于下限频率fl和上限频率fu之间。

三分之一倍频程谱是按逐级式频率进行分析的，它是由多个带通滤波器并联组成，为的是使这些带通滤波器的带宽覆盖整个分析频带。根据国籍电工委员会（IEC）的推荐，三分之一倍频程的中心频率为：

fc = 1000*10^(3n/30)Hz (n = ...-3,-2,-1,0,1,2,3...)

但在实际应用中，通常采用中心频率是其近似值。按照我国现行标准规定，中心频率为1Hz，1.25Hz，1.6Hz，2Hz，2.5Hz，3.15Hz，4Hz，5Hz，6.3Hz，8Hz，10Hz，...。可以看出，每隔三个中心频率，频率值增加一倍。三分之一倍频程的上下限频率以及中心频率之间的关系为：

fu/fl = 2^(1/3);fc/fl = 2^(1/6);fu/fc = 2^(1/6)

三分之一倍频程带宽为

△f = fu - fl

对于三分之一倍频程谱可以通过两种处理方法得到。一种方法是在整个分析频率从范围，按照不同的中心频率从定义对采样信号进行带通滤波，然后，计算出滤波后数据的均方值或均方根值（有效值），这样，变得到对应每个中心频率的功率谱值或幅值谱值。由于但分之一倍频程谱的滤波带与中心频率的比值是不变的，这种处理方法称为恒定百分比带宽滤波法。另一种方法首先对采样信号进行快速傅立叶变换，计算出功率谱或幅值谱，然后用功率谱或幅值谱的数据，计算每一个中心频率带宽内数矩的平均值，这样便处理得到三分之一倍频程谱值。这种方法显然要比第一种方法处理效率要高的多。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 三分之一倍频程处理
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear              % Remove items from workspace, freeing up system memory
clc                % Clear Command Window
clf                % Clear current figure window
close all hidden   % removal draws only those lines that are not obscured
                   % by other objects in the field of view.
format long
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 加载数据
fun = @(t) sin(50*2*pi*t) +sin(20*2*pi*t) + randn(size(t));
sf = 300;                       % 采样频率
n = 1000;                      % 数据长度
t = 0:1/sf:(n-1)/sf;            % 建立离散时间列向量
x = fun(t);                     % 产生时间序列数据

% 定义三分之一倍频程的中心频率
f = [1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.300 8.00];
fc = [f,10*f,100*f,1000*f,10000*f];
oc6 = 2^(1/6);        % 中心频率与下限频率的比值
nc = length(fc);      % 取中心频率总的长度     
nfft = 2^nextpow2(n); % 大于并最接近n的2的幂次方长度
a = fft(x,nfft);      % FFT变换
yc = zeros(1,nc);
for j = 1:nc   
      fl = fc(j)/oc6;  % 下限频率   
      fu = fc(j)*oc6;      % 上限频率
      nl = round(fl*nfft/sf+1);  % 下限频率对应的序号
      nu = round(fu*nfft/sf+1);  % 上限频率对应的序号
      if fu > sf/2     % 如果上限频率大于折叠频率则循环中断
         m = j-1; break
      end
    % 以每个中心频率段为通带进行带通频域滤波 
      b = zeros(1,nfft);
      b(nl:nu) = a(nl:nu);
      b(nfft-nu+1:nfft-nl+1) = a(nfft-nu+1:nfft-nl+1);
      c = ifft(b,nfft);
    %计算对应每个中心频率段的有效值
      yc(j) = sqrt(var(real(b(1:n))));
end
%绘制输入时程曲线图形
subplot(2,1,1);plot(t,x);
xlabel('时间　(s)');ylabel('加速度 (g)');grid on;                     
%绘制三分之一倍频程有效值图形
subplot(2,1,2); plot(fc(1:m),yc(1:m));
xlabel('频率 (Hz)');ylabel('有效值'); grid on;