hms2006 发表于 2007-8-14 17:30

解释下这段程序?

某天在本论坛上下了个程序,当时也没注意是什么,好象是和小波相关的,但后来随便演示了一下,还真不明白,到底是什么作用,因此请求大侠解释下这段程序。其中a1和 b1是小波分解后的的 逼近系数和细节系数。
%%%%%%%%%%%%%%%hilbert
fs=1000;
Y2=fft(a1,512);
Pyy2=Y2.*conj(Y2)/512;
f2=1000*(0:256)/512;
y=hilbert(a1);
ydata=abs(y);
y=y-mean(y);
nfft=1000;
p=abs(fft(ydata,nfft));
y1=hilbert(d1);
ydata1=abs(y1);
y1=y1-mean(y1);
nfft=1000;
p1=abs(fft(ydata1,nfft));


figure
subplot(311)
plot((0:nfft/2-1)/nfft*fs,p(1:nfft/2));
subplot(312)
plot((0:nfft/2-1)/nfft*fs,p1(1:nfft/2));
subplot(313)
plot(f2,Pyy2(1:257));

zhlong 发表于 2007-8-14 20:24

回复 #1 hms2006 的帖子

楼主最好能给出完整的原文。

hms2006 发表于 2007-8-14 21:01

好的,我把前面的一部分补上,其实那部分只是一个小波分解和重构程序,很普通,没什么的。
我只是不明白后面,也就是我贴出来的,又是FFT又是HILBERT,一会又求平均值到底有什么作用,这样处理的结果给出的是什么物理量或是物理意义。
%含噪的三角波与正弦波的组合
%生成正弦信号
clc;
clear;
close all;
N=1000;
t=0:0.001:0.6;
x=sin(2*pi*50*t)+sin(2*pi*120*t);
y=x+2*randn(size(t));
figure;
subplot(211)
plot(t,y,'LineWidth',2);
subplot(212)
plot(1000*t(1:50),y(1:50));
%xlabel('sample number n');
%ylabel('amplitude A');
%wavelet 1D decompose
=wavedec(y,7,'db10');
% reconstruct 1-7 layer approach coefficient.
a7=wrcoef('a',c,l,'db10',7);
a6=wrcoef('a',c,l,'db10',6);
a5=wrcoef('a',c,l,'db10',5);
a4=wrcoef('a',c,l,'db10',4);
a3=wrcoef('a',c,l,'db10',3);
a2=wrcoef('a',c,l,'db10',2);
a1=wrcoef('a',c,l,'db10',1);

% show approach coefficient
figure
subplot(7,1,1);
plot(a7,'LineWidth',2);
ylabel('a7');
subplot(7,1,2);
plot(a6,'LineWidth',2);
ylabel('a6');
subplot(7,1,3);
plot(a5,'LineWidth',2);
ylabel('a5');
subplot(7,1,4);
plot(a4,'LineWidth',2);
ylabel('a4');
subplot(7,1,5);
plot(a3,'LineWidth',2);
ylabel('a3');
subplot(7,1,6);
plot(a2,'LineWidth',2);
ylabel('a2');
subplot(7,1,7);
plot(a1,'LineWidth',2);
ylabel('a1');
xlabel('sample sequence n');
% reconstruction 1-7 layer details coefficient
d7=wrcoef('d',c,l,'db10',7);
d6=wrcoef('d',c,l,'db10',6);
d5=wrcoef('d',c,l,'db10',5);
d4=wrcoef('d',c,l,'db10',4);
d3=wrcoef('d',c,l,'db10',3);
d2=wrcoef('d',c,l,'db10',2);
d1=wrcoef('d',c,l,'db10',1);
% show detail coefficient
figure
subplot(7,1,1);
plot(d7,'LineWidth',2);
ylabel('d7');
subplot(7,1,2);
plot(d7,'LineWidth',2);
ylabel('d6');
subplot(7,1,3);
plot(d5,'LineWidth',2);
ylabel('db5');
subplot(7,1,4);
plot(d5,'LineWidth',2);
ylabel('d4');
subplot(7,1,5);
plot(d3,'LineWidth',2);
ylabel('d3');
subplot(7,1,6);
plot(d2,'LineWidth',2);
ylabel('d2');
subplot(7,1,7);
plot(d1,'LineWidth',2);
ylabel('d1');
xlabel('sample squence n');

zhlong 发表于 2007-8-14 22:08

fs=1000;
Y2=fft(a1,512);
Pyy2=Y2.*conj(Y2)/512;
f2=1000*(0:256)/512;
y=hilbert(a1);
ydata=abs(y);
y=y-mean(y);    %这句似乎应该为 ydata=ydata-mean(ydata);
nfft=1000;
p=abs(fft(ydata,nfft));%整个这段红色代码是求a1的包络谱

y1=hilbert(d1);
ydata1=abs(y1);
y1=y1-mean(y1); %这句似乎应该为 ydata1=ydata1-mean(ydata1);
nfft=1000;
p1=abs(fft(ydata1,nfft));%这段是求d1的包络谱


figure
subplot(311)
plot((0:nfft/2-1)/nfft*fs,p(1:nfft/2));
subplot(312)
plot((0:nfft/2-1)/nfft*fs,p1(1:nfft/2));
subplot(313)
plot(f2,Pyy2(1:257));
页: [1]
查看完整版本: 解释下这段程序?