倒频谱

sf=1000;

nfft=1000;

x=0:1/sf:1;

y=5*cos(2*pi*5*x).^2+3*cos(2*pi*100*x)+randn(size(x)); t=0:1/sfnfft-1)/sf;

nn=1:nfft/4;

subplot(3,1,1)

z=rceps(y);

plot(t(nn),abs(z(nn)));

title('z=rceps(y)')

ylabel('幅值');

grid on;

subplot(3,1,2)

yy= real(ifft(log(abs(fft(y)))));

plot(t(nn),abs(yy(nn)));

title('real(ifft(log(abs(fft(x)))))')

xlabel('时间(s)');

ylabel('幅值');

grid on;

subplot(3,1,3)

w=hanning(nfft);

zw=psd(y,nfft,sf,w,nfft/2);

zzw= real(ifft(log(abs(zw))));

plot(t(nn),zzw(nn));

grid on

为了方便楼主处理，上传了sunday.txt。下载后把后缀改为.wav。读取数

据的程序为：

waveFile='sunday.wav';

[speech, sf, nbits]=wavread(waveFile);

index1=4000;

nfft=1000;

index2=index1+nfft-1;

y=speech(index1:index2);

附件: sunday ( 30.74 K )

#4

dr

ea m _s ail or 工程师

精华 0 积分 208 帖子 104 水位 208

技

术

分

财

富

0 根据‘songzy41’提供的资源，编辑M 文件如下：

clear

clc

close all hidden

format long

waveFile='sunday.wav';

[speech, sf, nbits]=wavread(waveFile);

index1=4000;

nfft=1000;

index2=index1+nfft-1;

y=speech(index1:index2); t=0:1/sfnfft-1)/sf; nn=1:nfft/4; subplot(3,1,1) z=rceps(y); %MATLAB 提供的倒频谱函数 plot(t(nn),abs(z(nn))); title('z=rceps(y)') ylabel('幅值');

grid on;

subplot(3,1,2)

%实倒谱定义一：是通过对时域信号的傅里叶变换的幅值求自然对数，然后再做傅里叶逆变换。

yy= real(ifft(log(abs(fft(y)))));

plot(t(nn),abs(yy(nn)));

title('real(ifft(log(abs(fft(x)))))')

xlabel('时间(s)');

ylabel('幅值');

grid on;

subplot(3,1,3)

%实倒谱定义二：是通过对时域信号的功率谱密度的幅值求自然对数，然后再做傅里叶逆变换。

w=hanning(nfft);

zw=psd(y,nfft,sf,w,nfft/2);

zzw= real(ifft(log(abs(zw)))); plot(t(nn),zzw(nn)); grid on 1.得到信号的倒频谱图如下，从第一、二图可以看出，两图谱图一致，因此，断定MATLAB 中的倒频谱函数采用的是定义一，即是通过对时域信号的傅里叶变换的幅值求自然对数，然后再做傅里叶逆变换。这样理解对否？ 2.是否从你给的信号中：基频成分1/0.005875＝172Hz ；1/0.01163＝86Hz 。

我构造了一个信号，M文件如下：

clc;

clear;

fs=5000;

nfft=1024;

n=1:nfft/4; %(1,2,3, (256)

f=0:fs/nfft:fs/2-fs/nfft; %频率向量－横坐标

i=0;

for k=1:6

for t=0:1/fs:1

i=i+1;

x(1,i)=4*sin(2*pi*25*k*t)*sin(2*pi*350*t)+2.5*sin(2*pi*350*t);

end

end

subplot(3,1,1)

y=abs(fft(x,nfft)); %FFT

plot(f(n),y(n))

xlabel('Frequency');

ylabel('|FFT)|')

grid on

subplot(3,1,2)

w=hanning(nfft);

z=psd(x,nfft,fs,w,nfft/2); %PSD

plot(f(n),abs(z(n)));

xlabel('Frequency(Hz)');

ylabel('|psd()|')

grid on

subplot(3,1,3)

t=1000*(0:1/fsnfft-1)/fs); %倒频谱的时间向量－横坐标ms zw=rceps(x); %rceps

plot(t(n),abs(zw(n)));

xlabel('ms')

ylabel('|Cepstrum|');

grid on;

得到如下图：图三中1000/40＝25Hz