logcos 发表于 2008-5-21 14:59

1/2倍频程 如何用Matlab求解

老师给的要求:
读出一个Wav噪声 得到其原波形图
对声音信号进行傅里叶变换 得到其频谱图
再运行傅里叶变换 得到倒谱
得到声音信号的1/2倍频程图

我的写的程序如下: 不会救1/2倍频程

=wavread('d:\wav.wav');
t=length(y);
subplot(3,1,1)
plot(y);
legend('波形图');
xlabel('时间(s)');
ylabel('幅度');
%---------------
x=fft(y.*hamming(t));
fm=5000*length(x)/fs;
f=(0:fm)*fs/length(x);
subplot(3,1,2);
plot(f,20*log10(abs(x(1:length(f))+eps)));
legend('频谱图');
xlabel('频率Hz');
ylabel('频率幅度db');

c=fft(log(abs(x)+eps));
ms1=fs/1000;
ms20=fs/50;
q=(ms1:ms20)/fs;
subplot(3,1,3);
plot(q,abs(c(ms1:ms20)));
legend('倒谱图');
xlabel('倒频');
ylabel('倒频谱幅度');

请问对不,再请高手指点一下,1/2倍频程怎么求,谢谢~

w89986581 发表于 2008-5-21 15:18

论坛上有1/3倍频程的matlab程序,你搜索一下。把其中的频率间隔设置改一下即可,依次乘上2^(1/4)。什么领域惯用1/2倍频程呢?

[ 本帖最后由 无水1324 于 2008-5-27 09:55 编辑 ]

logcos 发表于 2008-5-21 16:05

谢谢楼上的朋友
找了一下

=wavread('d:\wav.wav');

N=2;
Fc=Fs/100;
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;
plot(butter(N,));

是这样的么?好像不对吧~
不知道老师为什么要让我做1/2的。

logcos 发表于 2008-5-21 16:35

1/2倍频程有错吧。请帮忙看看



前三个图是这样的

w89986581 发表于 2008-5-21 16:55

呵呵,我指的不是这个程序,你需要验证一下是否正确,而且根据1/2倍频程要求,需要将
f1 = Fc/(2^(1/6));
f2 = Fc*(2^(1/6));
中1/6换成1/4

w89986581 发表于 2008-5-21 16:56

程序不是写完就over了,还需要验证。最简单的方法就是构造一个单频函数加白噪声,用你的程序去分析,然后根据结果判断是否准确。

logcos 发表于 2008-5-21 16:59

谢谢楼上的指点,
实在是没有时间,以前少有接触这方面的。
着急要。

楼上的朋友能指点一下,简单一点的实现么?
谢谢~~

w89986581 发表于 2008-5-21 17:21

搜索到现成的程序,把1/6都改成1/4就可以拉,不能再简单拉,呵呵。最后自己用白噪声测试一下。

logcos 发表于 2008-5-21 17:27

我就找到上面这一个。好像不对。
楼上朋友请帮提供一个吧,感激不尽~
:handshake
:handshake

就差1/2倍频倍了,不然老师不让我毕业了。:@(
Q:38774566
等 待中。。。。。。。。

songzy41 发表于 2008-5-21 20:05

我把1/3倍频程滤波的程序改了一下,成1/2倍频程。1/3倍频程滤波器的频带是有国际标准的,而1/2倍频程没有标准。程序为
x=load('JGS-UD.txt');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%预处理
L=length(x);
amean=sum(x)/L;
a=x-amean;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sf=8000; %采样频率
f=[1.000    1.414    2.000    2.828    4.000    5.657    8.000   11.314   16.000   22.627...   
32.000   45.255   64.000   90.510   128.000   181.019   256.000   362.039   512.000   724.077];
fc=;
oc4=2^(1/4);
nc=length(fc);
n=length(x);
nfft=2^nextpow2(n);
freqStep=sf/nfft;
a1= fft(a,nfft);
for j=1:nc
fl=fc(j)/oc4;
fu=fc(j)*oc4;
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)=a1(nl:nu);
b(nfft-nu+1:nfft-nl+1)=a1(nfft-nu+1:nfft-nl+1);
c=ifft(b,nfft);
yc(j)=sqrt(var(real(c(1:n))));
end
plot(fc(1:m),yc(1:m));
xlabel('频率(Hz)');
ylabel('RMS(dB)');
grid on;
对语音信号只要把
x=load('JGS-UD.txt');
sf=8000; %采样频率
改为
=wavread('d:\wav.wav');

songzy41 发表于 2008-5-22 07:31

上帖说到,1/2倍频程没有标准,1/3倍频程和倍频程有标准,所以考虑到倍频程的标准,1/2倍频程的频率设置成如下的数值可能会更好:
f=[ 1.000    1.414    2.000    2.828    4.000    5.657    8.000   11.314   16.000   22.627...   
31.500   44.548   63.000   89.095   125.000   176.777   250.000   353.553   500.000   707.107];

logcos 发表于 2008-5-22 19:37

谢谢楼上的朋友,努力学习……
:handshake :handshake

《Matlab在振动信号处理中的运用》
页: [1]
查看完整版本: 1/2倍频程 如何用Matlab求解