MATLAB实现倍频程

yi_lu1025

如何在MATLAB实现倍频程？

Eminem

倍频程是声学里人的可听频率范围内，将声音的频谱进行一定规则的集中，变成有限的几个频点对应的强度，这样描述比较起来容易，是一种公约的描述形式。

使用1/3倍频程主要是因为人耳对声音的感觉，其频率分辨能力不是单一频率，而是频带，而1/3倍频程曾经被认为是比较符合人耳特性的频带划分方法，不过现在心理声学里提出了Critical Band这么个频带划分方法，听说更符合人耳特性，但1/3倍频程仍在广泛使用。

分析频谱时，对于连续谱而言，分析某频率点上的声功率是没有意义的，因此有必要统计某一频带内的声功率。对于频带划分，倍频程和1/3倍频程是常用的划分方法之一，它们都是相对恒定带宽，例如1/3倍频程的带宽是中心频率的23%。

声学及振动测量仪器中的倍频程及1/3倍频程滤波主要是用于对噪声或振动进行频谱分析用的,它们是一种等百分比带宽滤波器,与人耳的频谱分析特性相似。在噪声测量中，使用1/3oct主要是将噪声的频率分布情况更直观的表示出来。便于今后的工作开展。

百分比=(2^(m/2)-2^(-m/2))*100%

其中m就是几倍频程，1/3倍频程m等于1/3。

先要知道1/3倍频程的划分方法，相关的书和国标都有公式和现成的数据表格，然后，将时间域的声信号fft变换到频率域，对定义的每个1/3倍频带的声压计算等效连续声压级。这就是1/3倍频程声压级。

1. function [g,f] = oct3spec(B,A,Fs,Fc,s,n);
   
2. % OCT3SPEC Plots a one-third-octave filter characteristics.
   
3. %    OCT3SPEC(B,A,Fs,Fc) plots the attenuation of the filter defined by
   
4. %    B and A at sampling frequency Fs. Fc is the center frequency of
   
5. %    the one-third-octave filter. The plot covers one decade on both sides
   
6. %    of Fc.
   
7. %
   
8. %    OCT3SPEC(B,A,Fs,Fc,'ANSI',N) superposes the ANSI Order-N analog
   
9. %    specification for comparison. Default is N = 3.
   
10. %
    
11. %    OCT3SPEC(B,A,Fs,Fc,'IEC',N) superposes the characteristics of the
    
12. %    IEC 61260 class N specification for comparison. Default is N = 1.
    
13. %
    
14. %    [G,F] = OCT3SPEC(B,A,Fs,Fc) returns two 512-point vectors with
    
15. %    the gain (in dB) in G and logarithmically spaced frequencies in F.
    
16. %    The plot can then be obtained by SEMILOGX(F,G)
    
17. %                              
    
18. %    See also OCT3DSGN, OCTSPEC, OCTDSGN.
    
19. % Author: Christophe Couvreur, Faculte Polytechnique de Mons ( Belgium)
    
20. 
    
21. % Last modification: Sept. 4, 1997, 11:00am.
    
22. % References:
    
23. %    [1] ANSI S1.1-1986 (ASA 65-1986): Specifications for
    
24. %        Octave-Band and Fractional-Octave-Band Analog and
    
25. %        Digital Filters, 1993.
    
26. %    [2] IEC 61260 (1995-08):  Electroacoustics -- Octave-Band and
    
27. %        Fractional-Octave-Band Filters, 1995.   
    
28. if (nargin < 4) | (nargin > 6)
    
29.   error('Invalide number of input arguments.');
    
30. end
    
31. ansi = 0;
    
32. iec = 0;
    
33. if nargin > 4
    
34.   if strcmp(lower(s),'ansi')
    
35.     ansi = 1;
    
36.     if nargin == 5
    
37.       n = 3;
    
38.     end
    
39.   elseif strcmp(lower(s),'cei') | strcmp(lower(s),'iec')
    
40.     iec = 1;
    
41.      if nargin == 5
    
42.       n = 1
    
43.     end
    
44.     if (n < 0) | (n > 3)
    
45.       error('IEC class must be 0, 1, or 2');
    
46.     end
    
47.   end
    
48. end
    
49. N = 512;
    
50. pi = 3.14159265358979;
    
51. F = logspace(log10(Fc/10),log10(min(Fc*10,Fs/2)),N);
    
52. H = freqz(B,A,2*pi*F/Fs);
    
53. G = 20*log10(abs(H));
    
54. % Set output variables
    
55. if nargout ~= 0
    
56.   g = G; f = F;
    
57.   return
    
58. end
    
59. % Generate the plot
    
60. if (ansi)                       % ANSI Order-n specification
    
61.   f = logspace(log10(Fc/10),log10(Fc*10),N);
    
62.   f1 = Fc/(2^(1/6));
    
63.   f2 = Fc*(2^(1/6));
    
64.   Qr = Fc/(f2-f1);
    
65.   Qd = (pi/2/n)/(sin(pi/2/n))*Qr;
    
66.   Af = 10*log10(1+Qd^(2*n)*((f/Fc)-(Fc./f)).^(2*n));
    
67.   semilogx(F,G,f,-Af,'--');
    
68.   legend('Filter',['ANSI order-' int2str(n)],0);
    
69. elseif (iec)                                  % CEI specification
    
70.   semilogx(F,G);
    
71.   hold on
    
72.   if n == 0
    
73.     tolup =  [ .15 .15 .15 .15 .15 -2.3 -18.0 -42.5 -62 -75 -75 ];
    
74.     tollow = [ -.15 -.2 -.4 -1.1 -4.5 -realmax -inf -inf -inf -inf -inf ];
    
75.   elseif n == 1
    
76.     tolup =  [ .3 .3 .3 .3 .3 -2 -17.5 -42 -61 -70 -70 ];
    
77.     tollow = [ -.3 -.4 -.6 -1.3 -5 -realmax -inf -inf -inf -inf -inf ];
    
78.   elseif n == 2
    
79.     tolup =  [ .5 .5 .5 .5 .5 -1.6 -16.5 -41 -55 -60 -60  ];
    
80.     tollow = [ -.5 -.6 -.8 -1.6 -5.5 -realmax -inf -inf -inf -inf -inf ];
    
81.   end
    
82.   % Reference frequencies in base 2 system
    
83.   f = Fc * [1 1.02676 1.05594 1.08776 1.12246 1.12246 1.29565 1.88695 ...
    
84.          3.06955 5.43474 NaN ];  
    
85.   f(length(f)) = realmax;
    
86.   ff = Fc * [1 0.97394 0.94702 0.91932 0.89090 0.89090 0.77181 0.52996 ...
    
87.          0.32578 0.18400 NaN ];  
    
88.   ff(length(ff)) = realmin;
    
89.   semilogx(F,G,f,tolup,'--');
    
90.   semilogx(F,G,f,tollow,'--');
    
91.   semilogx(F,G,ff,tolup,'--');
    
92.   semilogx(F,G,ff,tollow,'--');
    
93.   hold off
    
94.   legend('Filter',['IEC class ' int2str(n)],0);
    
95. else
    
96.   semilogx(F,G);
    
97. end
    
98. xlabel('Frequency [Hz]'); ylabel('Gain [dB]');
    
99. title(['One-third-octave filter: Fc =',int2str(Fc),' Hz, Fs = ',int2str(Fs),' Hz']);
    
100. axis([Fc/10 Fc*10 -80 5]);
     
101. grid on
     
102. 
     
103. function [B,A] = oct3dsgn(Fc,Fs,N);
     
104. % OCT3DSGN  Design of a one-third-octave filter.
     
105. %    [B,A] = OCT3DSGN(Fc,Fs,N) designs a digital 1/3-octave filter with
     
106. %    center frequency Fc for sampling frequency Fs.
     
107. %    The filter is designed according to the Order-N specification
     
108. %    of the ANSI S1.1-1986 standard. Default value for N is 3.
     
109. %    Warning: for meaningful design results, center frequency used
     
110. %    should preferably be in range Fs/200 < Fc < Fs/5.
     
111. %    Usage of the filter: Y = FILTER(B,A,X).
     
112. %
     
113. %    Requires the Signal Processing Toolbox.
     
114. %
     
115. %    See also OCT3SPEC, OCTDSGN, OCTSPEC.
     
116. % Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
     
117. 
     
118. % Last modification: Aug. 25, 1997, 2:00pm.
     
119. % References:
     
120. %    [1] ANSI S1.1-1986 (ASA 65-1986): Specifications for
     
121. %        Octave-Band and Fractional-Octave-Band Analog and
     
122. %        Digital Filters, 1993.
     
123. if (nargin > 3) | (nargin < 2)
     
124.   error('Invalide number of arguments.');
     
125. end
     
126. if (nargin == 2)
     
127.   N = 3;
     
128. end
     
129. if (Fc > 0.88*(Fs/2))
     
130.   error('Design not possible. Check frequencies.');
     
131. end
     
132. % Design Butterworth 2Nth-order one-third-octave filter
     
133. % Note: BUTTER is based on a bilinear transformation, as suggested in [1].
     
134. pi = 3.14159265358979;
     
135. f1 = Fc/(2^(1/6));
     
136. f2 = Fc*(2^(1/6));
     
137. Qr = Fc/(f2-f1);
     
138. Qd = (pi/2/N)/(sin(pi/2/N))*Qr;
     
139. alpha = (1 + sqrt(1+4*Qd^2))/2/Qd;
     
140. W1 = Fc/(Fs/2)/alpha;
     
141. W2 = Fc/(Fs/2)*alpha;
     
142. [B,A] = butter(N,[W1,W2]);
     
143. 
     
144. function [p,f] = oct3bank(x);
     
145. % OCT3BANK Simple one-third-octave filter bank.
     
146. %    OCT3BANK(X) plots one-third-octave power spectra of signal vector X.
     
147. %    Implementation based on ANSI S1.11-1986 Order-3 filters.
     
148. %    Sampling frequency Fs = 44100 Hz. Restricted one-third-octave-band
     
149. %    range (from 100 Hz to 5000 Hz). RMS power is computed in each band
     
150. %    and expressed in dB with 1 as reference level.
     
151. %
     
152. %    [P,F] = OCT3BANK(X) returns two length-18 row-vectors with
     
153. %    the RMS power (in dB) in P and the corresponding preferred labeling
     
154. %    frequencies (ANSI S1.6-1984) in F.
     
155. %                              
     
156. %    See also OCT3DSGN, OCT3SPEC, OCTDSGN, OCTSPEC.
     
157. % Author: Christophe Couvreur, Faculte Polytechnique
     
158. 
     

复制代码

怪咖先生

matlab进行内插----就是在一组数种插入0
比如以前数组为A=[1 2 3 4 5]
现在可以内插后就编程A=[1 0 0 2 0 0 3 0 0 4 0 0 5 0 0],如果这个是时域的变化，时域扩大3倍，根据傅里叶变换性质就知道，频域缩小成原来的1/3，幅度也是原来的1/3，内插补零以后，经过fft，就可以达到1/3倍频，只需要再乘以系数3，就可以了
抽取---就是隔几个抽取一次
在频域中进行抽取也可以。
比如A=[1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 9
只需要编写n=1:(size(1)/3) B(n)=A(3*n),直接在频域进行抽取可以满足LZ的要求，B=[1 4 2 5 3 9]

Zoyie

好东西啊，谢谢各位大神~~~~

truleeee

怪咖先生 发表于 2016-10-10 09:06
matlab进行内插----就是在一组数种插入0
比如以前数组为A=[1 2 3 4 5]
现在可以内插后就编程A=[1 0 0 2 0 ...

这就是传说中的补零吗

yi_lu1025

Eminem 发表于 2016-10-8 13:34
倍频程是声学里人的可听频率范围内，将声音的频谱进行一定规则的集中，变成有限的几个频点对应的强度，这样 ...

谢谢

truleeee

Eminem 发表于 2016-10-8 13:34
倍频程是声学里人的可听频率范围内，将声音的频谱进行一定规则的集中，变成有限的几个频点对应的强度，这样 ...

试试看看好使不  感谢分享