信号处理方法 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

yangzj 发表于 2007-5-26 20:30

回复 #16 zhlong 的帖子

仔细看了下，确实有这个问题，资料中对这个功率谱的得来说的不清楚。

这样看起来的话是不是做短时分析后再取平均值。

zhlong 发表于 2007-5-26 20:34

回复 #17 yangzj 的帖子

^_^，看样子只有写那个资料的人最清楚了。
平均一下不知道表达的是什么意思？和直接FFT的区别呢？

lovelydeath 发表于 2007-5-26 20:35

我也是怎么都不懂资料中的方法
这是我原来的程序
function = sheji(file)
data=wavread(file);%语音文件读取
dat(1,:)=data(:,1);%取单路数据
ave=mean(dat);%计算直流分量
dat=dat-ave;%减去直流分量
w=hamming(256);%hamming 窗函数
dat=conv(dat,w);%卷积滤波
n=length(dat);%求dat长度
X=fft(dat);%进行FFT变换
%subplot(2,1,1);
%plot(abs(X));%画图
Y(:,floor(1:n/5))=X(:,floor(1:n/5));%取能量较集中的低频段
%subplot(2,1,2);
Y=abs(Y);
Y=Y.^2;%Y矩阵元素平方
plot(Y);%画图
这个明显我把信号全部都当成了平稳的做
而且我得出的东西就是尖尖的图我怀疑这个图没有任何意义

lovelydeath 发表于 2007-5-26 20:36

有没有可能是用PSD画的要怎么画呢

yangzj 发表于 2007-5-26 21:10

从这个资料附录的程序里找到他用的是
=spectrum(x,512,128,hamming(216),22050);

spectrum是老版本的用法，我还没找到它的用法，呵呵

zhlong 发表于 2007-5-26 21:18

你从哪里看到的？spectrum不是你那种用法啊？

help spectrum
SPECTRUM Spectral Estimation.
H = SPECTRUM.<ESTIMATOR> returns a spectral estimator in H of type
specified by ESTIMATOR.

There are three types of estimators: Power Spectral Density (PSD),
Mean-square Spectrum (MSS), and Pseudo Spectrum.Valid values for
<ESTIMATOR> for each of these three types are listed below. Type
"help spectrum.<ESTIMATOR>" to get help on a specific spectrum
estimator -e.g., "help spectrum.welch"):

----------------------------
Power Spectral Density (PSD)
----------------------------
burg - Burg                         periodogram - Periodogram
cov- Covariance                   welch    - Welch
mcov - Modified Covariance          yulear    - Yule-Walker AR
mtm- Thomson multitaper method (MTM)

To calculate the PSD you must create one of the estimators above and
pass it to one of the functions below along with the data.

Fs = 1000; t = 0:1/Fs:.296;
   x = cos(2*pi*t*200)+randn(size(t));
   h = spectrum.welch;                % Create a Welch spectral estimator.

   Hpsd = psd(h,x,'Fs',Fs);          % Calculate the PSD
   plot(Hpsd)                         % Plot the PSD.

yangzj 发表于 2007-5-26 21:25

我在网上搜了楼主给的资料，里面有附程序，那里spectrum用的是老版用法

function y=psdfun(x)%函数用来求信号的特征频率向量
=spectrum(x,512,128,hamming(216),22050);
m=0;
for k=2:length(p)-1;
   if p(k-1)<p(k)&p(k)>p(k+1)
   m=m+1;
   y(m)=f(k);
   end
end

function frez=eigfrez(x,N)%函数从某人的一组语音信号的特征频率向量中提取出标准特征向量frez，x为这组特征频率向量的细胞矩阵，N为信号数。
bin=zeros(N,60);n=zeros(N,60);y=zeros(N,60);
for i=1:N;
   for j=1:length(x{i});
   temp=fix(x{i}(1,j)/100);
   switch temp
         case 0
            bin(i,1)=1;n(i,1)=n(i,1)+1;y(i,1)=[(n(i,1)-1)*y(i,1)+x{i}(1,j)]/n(i,1);
         case 1
            bin(i,2)=1;n(i,2)=n(i,2)+1;y(i,2)=[(n(i,2)-1)*y(i,2)+x{i}(1,j)]/n(i,2);
         case 2
            bin(i,3)=1;n(i,3)=n(i,3)+1;y(i,3)=[(n(i,3)-1)*y(i,3)+x{i}(1,j)]/n(i,3);
         case 3
            bin(i,4)=1;n(i,4)=n(i,4)+1;y(i,4)=[(n(i,4)-1)*y(i,4)+x{i}(1,j)]/n(i,4);
         case 4
            bin(i,5)=1;n(i,5)=n(i,5)+1;y(i,5)=[(n(i,5)-1)*y(i,5)+x{i}(1,j)]/n(i,5);
         case 5
            bin(i,6)=1;n(i,6)=n(i,6)+1;y(i,6)=[(n(i,6)-1)*y(i,6)+x{i}(1,j)]/n(i,6);
         case 6
            bin(i,7)=1;n(i,7)=n(i,7)+1;y(i,7)=[(n(i,7)-1)*y(i,7)+x{i}(1,j)]/n(i,7);
         case 7
            bin(i,8)=1;n(i,8)=n(i,8)+1;y(i,8)=[(n(i,8)-1)*y(i,8)+x{i}(1,j)]/n(i,8);
         case 8
            bin(i,9)=1;n(i,9)=n(i,9)+1;y(i,9)=[(n(i,9)-1)*y(i,9)+x{i}(1,j)]/n(i,9);
         case 9
            bin(i,10)=1;n(i,10)=n(i,10)+1;y(i,10)=[(n(i,10)-1)*y(i,10)+x{i}(1,j)]/n(i,10);
         case 10
            bin(i,11)=1;n(i,11)=n(i,11)+1;y(i,11)=[(n(i,11)-1)*y(i,11)+x{i}(1,j)]/n(i,11);
      case 11
            bin(i,12)=1;n(i,12)=n(i,12)+1;y(i,12)=[(n(i,12)-1)*y(i,12)+x{i}(1,j)]/n(i,12);
      case 12
            bin(i,13)=1;n(i,13)=n(i,13)+1;y(i,13)=[(n(i,13)-1)*y(i,13)+x{i}(1,j)]/n(i,13);
         case 13
            bin(i,14)=1;n(i,14)=n(i,14)+1;y(i,14)=[(n(i,14)-1)*y(i,14)+x{i}(1,j)]/n(i,14);
         case 14
            bin(i,15)=1;n(i,15)=n(i,15)+1;y(i,15)=[(n(i,15)-1)*y(i,15)+x{i}(1,j)]/n(i,15);
         case 15
            bin(i,16)=1;n(i,16)=n(i,16)+1;y(i,16)=[(n(i,16)-1)*y(i,16)+x{i}(1,j)]/n(i,16);
         case 16
            bin(i,17)=1;n(i,17)=n(i,17)+1;y(i,17)=[(n(i,17)-1)*y(i,17)+x{i}(1,j)]/n(i,17);
         case 17
            bin(i,18)=1;n(i,18)=n(i,18)+1;y(i,18)=[(n(i,18)-1)*y(i,18)+x{i}(1,j)]/n(i,18);
         case 18
            bin(i,19)=1;n(i,19)=n(i,19)+1;y(i,19)=[(n(i,19)-1)*y(i,19)+x{i}(1,j)]/n(i,19);
         case 19
            bin(i,20)=1;n(i,20)=n(i,20)+1;y(i,20)=[(n(i,20)-1)*y(i,20)+x{i}(1,j)]/n(i,20);

      case 20
            bin(i,21)=1;n(i,21)=n(i,21)+1;y(i,21)=[(n(i,21)-1)*y(i,21)+x{i}(1,j)]/n(i,21);
      case 21
            bin(i,22)=1;n(i,22)=n(i,22)+1;y(i,22)=[(n(i,22)-1)*y(i,22)+x{i}(1,j)]/n(i,22);
      case 22
            bin(i,23)=1;n(i,23)=n(i,23)+1;y(i,23)=[(n(i,23)-1)*y(i,23)+x{i}(1,j)]/n(i,23);
         case 23
            bin(i,24)=1;n(i,24)=n(i,24)+1;y(i,24)=[(n(i,24)-1)*y(i,24)+x{i}(1,j)]/n(i,24);
         case 24
            bin(i,25)=1;n(i,25)=n(i,25)+1;y(i,25)=[(n(i,25)-1)*y(i,25)+x{i}(1,j)]/n(i,25);
         case 25
            bin(i,26)=1;n(i,26)=n(i,26)+1;y(i,26)=[(n(i,26)-1)*y(i,26)+x{i}(1,j)]/n(i,26);
         case 26
            bin(i,27)=1;n(i,27)=n(i,27)+1;y(i,27)=[(n(i,27)-1)*y(i,27)+x{i}(1,j)]/n(i,27);
         case 27
            bin(i,28)=1;n(i,28)=n(i,28)+1;y(i,28)=[(n(i,28)-1)*y(i,28)+x{i}(1,j)]/n(i,28);
         case 28
            bin(i,29)=1;n(i,29)=n(i,29)+1;y(i,29)=[(n(i,29)-1)*y(i,29)+x{i}(1,j)]/n(i,29);
         case 29
            bin(i,30)=1;n(i,30)=n(i,30)+1;y(i,30)=[(n(i,30)-1)*y(i,30)+x{i}(1,j)]/n(i,30);
      case 30
            bin(i,31)=1;n(i,31)=n(i,31)+1;y(i,31)=[(n(i,31)-1)*y(i,31)+x{i}(1,j)]/n(i,31);
      case 31
            bin(i,32)=1;n(i,32)=n(i,32)+1;y(i,32)=[(n(i,32)-1)*y(i,32)+x{i}(1,j)]/n(i,32);
      case 32
            bin(i,33)=1;n(i,33)=n(i,33)+1;y(i,33)=[(n(i,33)-1)*y(i,33)+x{i}(1,j)]/n(i,33);
         case 33
            bin(i,34)=1;n(i,34)=n(i,34)+1;y(i,34)=[(n(i,34)-1)*y(i,34)+x{i}(1,j)]/n(i,34);
         case 34
            bin(i,35)=1;n(i,35)=n(i,35)+1;y(i,35)=[(n(i,35)-1)*y(i,35)+x{i}(1,j)]/n(i,35);
         case 35
            bin(i,36)=1;n(i,36)=n(i,36)+1;y(i,36)=[(n(i,36)-1)*y(i,36)+x{i}(1,j)]/n(i,36);
         case 36
            bin(i,37)=1;n(i,37)=n(i,37)+1;y(i,37)=[(n(i,37)-1)*y(i,37)+x{i}(1,j)]/n(i,37);
         case 37
            bin(i,38)=1;n(i,38)=n(i,38)+1;y(i,38)=[(n(i,38)-1)*y(i,38)+x{i}(1,j)]/n(i,38);
         case 38
            bin(i,39)=1;n(i,39)=n(i,39)+1;y(i,39)=[(n(i,39)-1)*y(i,39)+x{i}(1,j)]/n(i,39);
         case 39
            bin(i,40)=1;n(i,40)=n(i,40)+1;y(i,40)=[(n(i,40)-1)*y(i,40)+x{i}(1,j)]/n(i,40);
         case40
            bin(i,41)=1;n(i,41)=n(i,41)+1;y(i,41)=[(n(i,41)-1)*y(i,41)+x{i}(1,j)]/n(i,41);
         case 41
            bin(i,42)=1;n(i,42)=n(i,42)+1;y(i,42)=[(n(i,42)-1)*y(i,42)+x{i}(1,j)]/n(i,42);
         case 42
            bin(i,43)=1;n(i,43)=n(i,43)+1;y(i,43)=[(n(i,43)-1)*y(i,43)+x{i}(1,j)]/n(i,43);
         case 43
            bin(i,44)=1;n(i,44)=n(i,44)+1;y(i,44)=[(n(i,44)-1)*y(i,44)+x{i}(1,j)]/n(i,44);
         case 44
            bin(i,45)=1;n(i,45)=n(i,45)+1;y(i,45)=[(n(i,45)-1)*y(i,45)+x{i}(1,j)]/n(i,45);
         case 45
            bin(i,46)=1;n(i,46)=n(i,46)+1;y(i,46)=[(n(i,46)-1)*y(i,46)+x{i}(1,j)]/n(i,46);
         case 46
            bin(i,47)=1;n(i,47)=n(i,47)+1;y(i,47)=[(n(i,47)-1)*y(i,47)+x{i}(1,j)]/n(i,47);
         case 47
            bin(i,48)=1;n(i,48)=n(i,48)+1;y(i,48)=[(n(i,48)-1)*y(i,48)+x{i}(1,j)]/n(i,48);
         case 48
            bin(i,49)=1;n(i,49)=n(i,49)+1;y(i,49)=[(n(i,49)-1)*y(i,49)+x{i}(1,j)]/n(i,49);
         case 49
            bin(i,50)=1;n(i,50)=n(i,50)+1;y(i,50)=[(n(i,50)-1)*y(i,50)+x{i}(1,j)]/n(i,50);
         case 50
            bin(i,51)=1;n(i,51)=n(i,51)+1;y(i,51)=[(n(i,51)-1)*y(i,51)+x{i}(1,j)]/n(i,51);
      case 51
            bin(i,52)=1;n(i,52)=n(i,52)+1;y(i,52)=[(n(i,52)-1)*y(i,52)+x{i}(1,j)]/n(i,52);
      case 52
            bin(i,53)=1;n(i,53)=n(i,53)+1;y(i,53)=[(n(i,53)-1)*y(i,53)+x{i}(1,j)]/n(i,53);
         case 53
            bin(i,54)=1;n(i,54)=n(i,54)+1;y(i,54)=[(n(i,54)-1)*y(i,54)+x{i}(1,j)]/n(i,54);
         case 54
            bin(i,55)=1;n(i,55)=n(i,55)+1;y(i,55)=[(n(i,55)-1)*y(i,55)+x{i}(1,j)]/n(i,55);
         case 55
            bin(i,56)=1;n(i,56)=n(i,56)+1;y(i,56)=[(n(i,56)-1)*y(i,56)+x{i}(1,j)]/n(i,56);
         case 56
            bin(i,57)=1;n(i,57)=n(i,57)+1;y(i,57)=[(n(i,57)-1)*y(i,57)+x{i}(1,j)]/n(i,57);
         case 57
            bin(i,58)=1;n(i,58)=n(i,58)+1;y(i,58)=[(n(i,58)-1)*y(i,58)+x{i}(1,j)]/n(i,58);
         case 58
            bin(i,59)=1;n(i,59)=n(i,59)+1;y(i,59)=[(n(i,59)-1)*y(i,59)+x{i}(1,j)]/n(i,59);
         case 59
            bin(i,60)=1;n(i,60)=n(i,60)+1;y(i,60)=[(n(i,60)-1)*y(i,60)+x{i}(1,j)]/n(i,60);
   end
end
end
t=sum(bin);m=0;
for k=1:60;
if t(1,k)>0;
      m=m+1;
      frez(m)=sum(y(:,k))/t(1,k);
end
end
function d=distance(from,to)%函数用来求特征向量from与标准向量to的距离d
Nf=length(from);
Nt=length(to);
for i=1:Nf;
iffrom(i)<6000;
   j=j+1;
   dis(j)=min((to-from(i)*ones(1,Nt)).^2);
end
end
d=sum(dis)/(j*10^4);

zhlong 发表于 2007-5-26 21:31

辛苦yangzj版主了

yangzj 发表于 2007-5-26 21:47

没找到这个函数的老版用法的说明。
调试运行了下，=spectrum(x,nfft,noverlap,window,fs);
各参数的含义：x为分析序列；
                  nfft为FFT点数；
                  noverlap是重叠采样点数；
                  window为窗，长度可自定
返回功率谱
做法还是重叠分段，分别求功率谱，再做平均

[ 本帖最后由 yangzj 于 2007-5-26 21:49 编辑 ]

lovelydeath 发表于 2007-5-26 22:12

是不是直接运行上面的就可以了？
我在台湾网上找到了这个函数FFTONESIDE
这是他的使用例子
% 顯示一個語音音框的單邊頻譜
=wavread('welcome.wav');
signal=y(2047:2047+237-1);
=fftOneSide(signal, fs, 1);
斑竹大人看看啊这个得出的中间第二个图可以不？
M函数见附件

[ 本帖最后由 lovelydeath 于 2007-5-26 22:16 编辑 ]

yangzj 发表于 2007-5-26 22:25

回复 #26 lovelydeath 的帖子

这个只是做一帧的单边谱，如果按你上面文章里程序的做法的话，他是做多帧（帧间可有重叠）再做平均

yangzj 发表于 2007-5-26 22:36

你可参照我上面对spectrum的说明来做：

令index=1:NWindow;设信号为x;

1、令y＝x(index),并对y加长度为NWindow的窗；
2、做yf=fft(y,NFFT);
3、得到y的单边功率谱: yP=(abs(yf(1:NFFT/2))).^2/(NFFT);
4、令index=index+NWindow－Noverlap（Noverlap为重叠采样点数）
重复1－4，直到到x终点；

最后把所有的yP做平均。

大概是这样，具体你实现下

lovelydeath 发表于 2007-5-26 22:41

回复 #27 yangzj 的帖子

我就用那个FFTONESIDE 提取出现波峰的频率做特征模版可不可以呢？

zhlong 发表于 2007-5-26 22:55

回复 #29 lovelydeath 的帖子

那个函数其实就是Fourier频谱分析吧，你直接用的话好像就跟你前面说的短时傅立叶分析没有关系了。

lovelydeath 发表于 2007-5-27 00:27

谢谢大家帮忙 :loveliness:

页: 1 [2] 3

声振论坛's Archiver

回复 #16 zhlong 的帖子

回复 #17 yangzj 的帖子

回复 #26 lovelydeath 的帖子

回复 #27 yangzj 的帖子

回复 #29 lovelydeath 的帖子