clc;
clear;
close all;
% 当前延拓模式是补零
oldmode = dwtmode('zpd');
%% 小波分解与重构
%装载信号
data = xlsread('1956-2011张家山月runoff.xlsx','sheet1');
s = data(12,:);
% index = 1:672;
% ts = 1;
% fs = 1/ts;
% t=1:ts:672;
% N = length(s);
% %画出原始信号
% figure(1);
% plot(index,s);
% ylable = ('jingliu');
% xlable = ('xuhao n');
%多尺度一维分解
[c,l] = wavedec(s,5,'db5');
% 重构第1～5层逼近系数
% a5 = wrcoef('a',c,l,'db5',5);
% a6 = wrcoef('a',c,l,'db5',6);
a5 = wrcoef('a',c,l,'db5',5);
a4 = wrcoef('a',c,l,'db5',4);
a3 = wrcoef('a',c,l,'db5',3);
a2 = wrcoef('a',c,l,'db5',2);
a1 = wrcoef('a',c,l,'db5',1);
% % 显示逼近系数
% figure
% subplot(6,1,1);plot(a6);ylabel('a6');
% subplot(6,1,2);plot(a5);ylabel('a5');
% subplot(6,1,3);plot(a4);ylabel('a4');
% subplot(6,1,4);plot(a3);ylabel('a3');
% subplot(6,1,5);plot(a2);ylabel('a2');
% subplot(6,1,6);plot(a1);ylabel('a1');
% xlabel('时间 t/s');
figure
subplot(5,1,1);plot(a5);ylabel('a5');
subplot(5,1,2);plot(a4);ylabel('a4');
subplot(5,1,3);plot(a3);ylabel('a3');
subplot(5,1,4);plot(a2);ylabel('a2');
subplot(5,1,5);plot(a1);ylabel('a1');
xlabel('时间 t/s');
% 重构第1～5层细节系数
% d6 = wrcoef('a',c,l,'db5',6);
d5 = wrcoef('d',c,l,'db5',5);
d4 = wrcoef('d',c,l,'db5',4);
d3 = wrcoef('d',c,l,'db5',3);
d2 = wrcoef('d',c,l,'db5',2);
d1 = wrcoef('d',c,l,'db5',1);
% 显示细节系数
figure
subplot(5,1,1);plot(d5);ylabel('d5');
subplot(5,1,2);plot(d4);ylabel('d4');
subplot(5,1,3);plot(d3);ylabel('d3');
subplot(5,1,4);plot(d2);ylabel('d2');
subplot(5,1,5);plot(d1);ylabel('d1');
xlabel('时间 t/年');
% figure
% subplot(6,1,1);plot(d6);ylabel('d6');
% subplot(6,1,2);plot(d5);ylabel('d5');
% subplot(6,1,3);plot(d4);ylabel('d4');
% subplot(6,1,4);plot(d3);ylabel('d3');
% subplot(6,1,5);plot(d2);ylabel('d2');
% subplot(6,1,6);plot(d1);ylabel('d1');
% xlabel('时间 t/年');
xxk = [];
xxk(:,1) = a1;
xxk(:,2) = a2;
xxk(:,3) = a3;
xxk(:,4) = a4;
xxk(:,5) = a5;
xxk(:,6) = d1;
xxk(:,7) = d2;
xxk(:,8) = d3;
xxk(:,9) = d4;
xxk(:,10) = d5;
% xxk(:,11) = d6;
% xxk(:,12) = a6;
xlswrite('xiaobo.xlsx',xxk);
% 第1层细节信号的包络谱
yh = hilbert(d3);
aabs = abs(yh); % 包络的绝对值
aabs = aabs - mean(aabs);
aangle = unwrap(angle(yh)); % 包络的相位
af = diff(aangle)/2/pi; % 包络的瞬时频率，差分代替微分计算
% NFFT = 2^nextpow2(N);
NFFT = 2^nextpow2(1024*4); % 改善栅栏效应
f = fs*linspace(0,1,NFFT);
YH = fft(yh, NFFT)/N; % Hilbert变换复信号的频谱
A = fft(aabs, NFFT)/N; % 包络的频谱
figure
plot(f,abs(YH))
title('原始复信号的Hilbert谱')
xlabel('频率f (Hz)')
ylabel('|YH(f)|')
figure
subplot(221)
plot(t, aabs')
title('包络的绝对值')
legend('包络分析结果', '真实包络')
subplot(222)
plot(t, aangle)
title('调制信号的相位')
subplot(223)
plot(t(1:end-1), af*fs)
title('调制信号的瞬时频率')
subplot(224)
plot(f,abs(A))
title('包络的频谱')
xlabel('频率f (Hz)')
ylabel('|A(f)|')
% 恢复延拓模式
dwtmode(oldmode);
% % % c = wavedec2(a,3,'db1');
% % %提取系数
% % cA3 = appcoef(C,L,'db1',3);
% % cD3 = detcoef(C,L,3);
% % cD2 = detcoef(C,L,2);
% % cD1 = detcoef(C,L,1);
% %
% % %重构系数
% % A3 = wrcoef('a',C,L,'db1',3);
% % D1 = wrcoef('d',C,L,'db1',2);
% % D2 = wrcoef('d',C,L,'db1',1);
% % D3 = wrcoef('d',C,L,'db1',2);
% %
% % %重构原始信号
% % A0 = waverec(C,L,'db1');
% %
% % %重构最大误差
% % Err = max(abs(a-A0));
% %
