matlab倍频程设计,[转载]1/3倍频程及Matlab程序实现

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

使用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倍频程声压级。

function [g,f] =

oct3spec(B,A,Fs,Fc,s,n);

% OCT3SPEC Plots a

one-third-octave filter characteristics.

%

OCT3SPEC(B,A,Fs,Fc) plots the attenuation of the filter defined

by

%

B and A at sampling frequency Fs. Fc is the center frequency

of

%

the one-third-octave filter. The plot covers one decade on both

sides

%

of Fc.

%

%

OCT3SPEC(B,A,Fs,Fc,'ANSI',N) superposes the ANSI Order-N

analog

%

specification for comparison. Default is N = 3.

%

%

OCT3SPEC(B,A,Fs,Fc,'IEC',N) superposes the characteristics of

the

%

IEC 61260 class N specification for comparison. Default is N =

1.

%

%

[G,F] = OCT3SPEC(B,A,Fs,Fc) returns two 512-point vectors

with

%

the gain (in dB) in G and logarithmically spaced frequencies in

F.

%

The plot can then be obtained by SEMILOGX(F,G)

%

%

See also OCT3DSGN, OCTSPEC, OCTDSGN.

% Author: Christophe

Couvreur, Faculte Polytechnique de Mons (

Belgium)

%

couvreur@thor.fpms.ac.be

% Last modification:

Sept. 4, 1997, 11:00am.

%

References:

%

[1] ANSI S1.1-1986 (ASA 65-1986): Specifications

for

%

Octave-Band and Fractional-Octave-Band Analog and

%

Digital Filters, 1993.

%

[2] IEC 61260 (1995-08): Electroacoustics --

Octave-Band and

%

Fractional-Octave-Band Filters, 1995.

if (nargin

< 4) | (nargin > 6)

error('Invalide

number of input arguments.');

end

ansi =

0;

iec =

0;

if nargin

> 4

if

strcmp(lower(s),'ansi')

ansi = 1;

if nargin == 5

n = 3;

end

elseif

strcmp(lower(s),'cei') | strcmp(lower(s),'iec')

iec = 1;

if nargin == 5

n = 1

end

if (n < 0) | (n >

3)

error('IEC class must be 0, 1, or 2');

end

end

end

N =

512;

pi =

3.14159265358979;

F =

logspace(log10(Fc/10),log10(min(Fc*10,Fs/2)),N);

H =

freqz(B,A,2*pi*F/Fs);

G =

20*log10(abs(H));

% Set output

variables

if nargout ~=

0

g = G; f =

F;

return

end

% Generate the

plot

if (ansi)

% ANSI Order-n specification

f =

logspace(log10(Fc/10),log10(Fc*10),N);

f1 =

Fc/(2^(1/6));

f2 =

Fc*(2^(1/6));

Qr =

Fc/(f2-f1);

Qd =

(pi/2/n)/(sin(pi/2/n))*Qr;

Af =

10*log10(1+Qd^(2*n)*((f/Fc)-(Fc./f)).^(2*n));

semilogx(F,G,f,-Af,'--');

legend('Filter',['ANSI order-' int2str(n)],0);

elseif (iec)

% CEI specification

semilogx(F,G);

hold

on

if n ==

0

tolup = [

.15 .15 .15 .15 .15 -2.3 -18.0 -42.5 -62 -75 -75

];

tollow = [ -.15 -.2 -.4 -1.1 -4.5 -realmax -inf -inf -inf -inf -inf

];

elseif n ==

1

tolup = [

.3 .3 .3 .3 .3 -2 -17.5 -42 -61 -70 -70 ];

tollow = [ -.3 -.4 -.6 -1.3 -5 -realmax -inf -inf -inf -inf -inf

];

elseif n ==

2

tolup = [

.5 .5 .5 .5 .5 -1.6 -16.5 -41 -55 -60 -60 ];

tollow = [ -.5 -.6 -.8 -1.6 -5.5 -realmax -inf -inf -inf -inf -inf

];

end

% Reference

frequencies in base 2 system

f = Fc * [1 1.02676

1.05594 1.08776 1.12246 1.12246 1.29565 1.88695

...

3.06955

5.43474 NaN ];

f(length(f)) =

realmax;

ff = Fc * [1 0.97394

0.94702 0.91932 0.89090 0.89090 0.77181 0.52996

...

0.32578

0.18400 NaN ];

ff(length(ff)) =

realmin;

semilogx(F,G,f,tolup,'--');

semilogx(F,G,f,tollow,'--');

semilogx(F,G,ff,tolup,'--');

semilogx(F,G,ff,tollow,'--');

hold

off

legend('Filter',['IEC

class ' int2str(n)],0);

else

semilogx(F,G);

end

xlabel('Frequency

[Hz]'); ylabel('Gain [dB]');

title(['One-third-octave filter: Fc

=',int2str(Fc),' Hz, Fs = ',int2str(Fs),' Hz']);

axis([Fc/10 Fc*10 -80

5]);

grid on

function [B,A] =

oct3dsgn(Fc,Fs,N);

%

OCT3DSGN

Design of a one-third-octave filter.

%

[B,A] = OCT3DSGN(Fc,Fs,N) designs a digital 1/3-octave filter

with

%

center frequency Fc for sampling frequency Fs.

%

The filter is designed according to the Order-N

specification

%

of the ANSI S1.1-1986 standard. Default value for N is

3.

%

Warning: for meaningful design results, center frequency

used

%

should preferably be in range Fs/200 < Fc

< Fs/5.

%

Usage of the filter: Y = FILTER(B,A,X).

%

%

Requires the Signal Processing Toolbox.

%

%

See also OCT3SPEC, OCTDSGN, OCTSPEC.

% Author: Christophe

Couvreur, Faculte Polytechnique de Mons (Belgium)

%

couvreur@thor.fpms.ac.be

% Last modification:

Aug. 25, 1997, 2:00pm.

%

References:

%

[1] ANSI S1.1-1986 (ASA 65-1986): Specifications

for

%

Octave-Band

and Fractional-Octave-Band Analog and

%

Digital Filters, 1993.

if (nargin

> 3) | (nargin < 2)

error('Invalide

number of arguments.');

end

if (nargin ==

2)

N =

3;

end

if (Fc

> 0.88*(Fs/2))

error('Design not

possible. Check frequencies.');

end

% Design Butterworth

2Nth-order one-third-octave filter

% Note: BUTTER is

based on a bilinear transformation, as suggested in

[1].

pi =

3.14159265358979;

f1 =

Fc/(2^(1/6));

f2 =

Fc*(2^(1/6));

Qr =

Fc/(f2-f1);

Qd =

(pi/2/N)/(sin(pi/2/N))*Qr;

alpha = (1 +

sqrt(1+4*Qd^2))/2/Qd;

W1 =

Fc/(Fs/2)/alpha;

W2 =

Fc/(Fs/2)*alpha;

[B,A] = butter(N,[W1,W2]);

function [p,f] =

oct3bank(x);

% OCT3BANK Simple

one-third-octave filter bank.

%

OCT3BANK(X) plots one-third-octave power spectra of signal vector

X.

%

Implementation based on ANSI S1.11-1986 Order-3

filters.

%

Sampling frequency Fs = 44100 Hz. Restricted

one-third-octave-band

%

range (from 100 Hz to 5000 Hz). RMS power is computed in each

band

%

and expressed in dB with 1 as reference level.

%

%

[P,F] = OCT3BANK(X) returns two length-18 row-vectors

with

%

the RMS power (in dB) in P and the corresponding preferred

labeling

%

frequencies (ANSI S1.6-1984) in F.

%

%

See also OCT3DSGN, OCT3SPEC, OCTDSGN, OCTSPEC.

% Author: Christophe

Couvreur, Faculte Polytechnique