📄 matlab_wavelet.txt
字号:
第 1 楼发表于 2006-3-31 08:34 资料 个人空间 短消息
[分享]个人收集的一些关于小波分析的matlab程序
都是从网上收集来的,由于时间比较久,处处都忘记了,如果是谁的原创请和我联系,
我在帖子中标出来的
内容比较多,将会逐步贴出来
提升法97经典程序 (二楼)
2代小波示意程序 (三楼)
二代小波漫谈 (四楼)
小波滤波器构造和消噪程序(2个) (五楼)
小波谱分析mallat算法经典程序 (六楼)
2维小波变换经典程序 (七楼)
基于LeventCodes平台的小波去噪程序包 (十一楼)
连续小波和离散小波分析的应用实例(十二楼)
小波插值与小波构造(3个程序)(十三楼)
采用多孔trous算法(undecimated wavelet transform)实现小波变换(十四楼)
Daubechies小波基的构造(十五楼)
消失矩作用的程序(二十三楼)
平移变换平移法(cycle_spinning)消除gibbs效应 (二十四楼)
第 2 楼发表于 2006-3-31 08:36 资料 个人空间 短消息
提升法97经典程序
[Copy to clipboard]CODE:
%% 本程序实现任意偶数大小图像第二代双正交97提升小波变换
%% 注1: 采用标准正交方法,对行列采用不同矩阵(和matlab里不同)
%% 注2: 为了保证正交,所有边界处理,全部采用循环处理
%% 注3: 正交性验证,将单位阵带入函数,输出仍是单位阵(matlab不具有此性质)
%% 注4: 此程序是矩阵实现,所以图像水平分量和垂直分量估计被交换位置
%% 注5: 此程序实现的是类小波(wavelet-like)变换,是介于小波包变换与小波变换之间的变换
%% 注6: 此程序每层变换相对原图像矩阵,产生的矩阵都是正交阵,这和小波包一致
%% 注7: 但小波变换每层产生的矩阵,是相对每个待分解子块的正交矩阵,而不是原图像的正交矩阵
%% 注8: 且小波变换产生的正交矩阵维数,随分解层数2分减少
%% 注9: 提升系数可以在MATLAB7.0以上版本,用liftwave('9.7')获取,这里直接给出,考虑兼容性
%% 注10:由于MATLAB数组下标从1开始,所以注意奇偶序列的变化
%% 注11:d为对偶上升,即预测;p为原上升,即更新 %% 编程人 沙威 安徽大学
%% 编程时间 2004年12月18日 %% x输入图像,y输出图像
%% flag_trans为正变换或反变换标志,0执行正变换,1执行反变换
%% flag_max,是否最大层数变换标志,0执行用户设定层数,1执行最大层数变换
%% layer,用户层数设置(小于最大层) function y=db97(x,flag_trans,flag_max,layer); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1.输入参数检查 % 矩阵维数判断
[sa,sb]=size(x); if (sa~=sb) % 防止非图像数据
errordlg('非图像数据!');
error('非图像数据!');
end; % 变换标志判断
[sa,sb]=size(flag_trans);
if ((sa~=1) | (sb~=1)) % 变换标志错误
errordlg('变换标志错误!');
error('变换标志错误!');
end; if ((flag_trans~=1) & (flag_trans~=0)) % 变换标志错误
errordlg('变换标志错误!');
error('变换标志错误!');
end; % 最大层数标志判断
[sa,sb]=size(flag_max);
if ((sa~=1) | (sb~=1)) % 最大层数标志错误
errordlg('最大层数标志错误!');
error('最大层数标志错误!');
end; if ((flag_max~=1) & (flag_max~=0)) % 最大层数标志错误
errordlg('最大层数标志错误!');
error('最大层数标志错误!');
end; % 用户设置层数判断
if (flag_max~=1) [sa,sb]=size(layer);
if ((sa~=1) | (sb~=1)) % 层数设置错误
errordlg('层数设置错误!');
error('层数设置错误!');
end; if (flag_max<0) % 层数设置错误
errordlg('层数设置错误!');
error('层数设置错误!');
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2.提升系数确定
% t1=liftwave('9.7'); % 获取提升系数(MATLAB7.0以后) d1=[-1.586100000000000e+000,-1.586134342069360e+000];
p1=[1.079600000000000e+000,-5.298011857188560e-002];
d2=[-8.829110755411875e-001,-8.829110755411875e-001];
p2=[4.435068520511142e-001,1.576123746148364e+000];
d3=-8.698644516247808e-001;
p3=-1.149604398860242e+000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.分解层数确定
% 采用用户输入和自动给出最大层数两种方法 N=length(x); % 矩阵大小
S=N; % 变量
s=log2(N); % 最大循环次数
n1=N/2; % 初始一半矩阵大小
n2=N; % 初始矩阵大小
u=0; % 初始值 % 对非2的整数幂大小图像确定最大分解层数
for ss=1:s
if (mod(S,2)==0)
u=u+1;
S=S/2;
end;
end;
u=u-1; % 分解最大层数减1(后面的边界处理造成) % 最大层数确定
if (flag_max==0) % 手动输入
T=layer; % 用户输入值
else % 自动确定最大层数
T=u; % 分解最大层数
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4.最大层数和图像大小检查 if (T>u) % 防止用户层数越界
errordlg('已超过最大分解层数!或者非偶数大小图像!');
error('已超过最大分解层数!或者非偶数大小图像!');
end; if (mod(N,2)~=0) % 防止图像大小错误
errordlg('非偶数大小图像!');
error('非偶数大小图像!');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5.提升法正变换 if (flag_trans==0)
for time=1:T; % 行正变换
% d;
x1(n1,:)=x(n2,:)+d1(2)*x(n2-1,:)+d1(1)*x(1,:);
x1([1:n1-1],:)=x([2:2:n2-2],:)+d1(2)*x([1:2:n2-3],:)+d1(1)*x([3:2:n2-1],:);
% p;
x(1,:)=x(1,:)+p1(2)*x1(n1,:)+p1(1)*x1(1,:);
x([2:n1],:)=x([3:2:n2-1],:)+p1(2)*x1([1:n1-1],:)+p1(1)*x1([2:n1],:);
x([n1+1:n2],:)=x1([1:n1],:);
% d;
x(n1+1,:)=x(n1+1,:)+d2(2)*x(n1,:)+d2(1)*x(1,:);
x([n1+2:n2],:)=x([n1+2:n2],:)+d2(2)*x([1:n1-1],:)+d2(1)*x([2:n1],:);
% p;
x(n1,:)=x(n1,:)+p2(2)*x(n1+1,:)+p2(1)*x(n1+2,:);
x(n1-1,:)=x(n1-1,:)+p2(2)*x(n2,:)+p2(1)*x(n1+1,:);
x([1:n1-2],:)=x([1:n1-2],:)+p2(2)*x([n1+2:n2-1],:)+p2(1)*x([n1+3:n2],:);
% 归一
x([1:n1],:)=p3*x([1:n1],:);
x([n1+1:n2],:)=d3*x([n1+1:n2],:); clear x1;
% 列正变换
% d;
x1(:,[1:n1])=x(:,[2:2:n2]);
% p;
x(:,1)=x(:,1)-d1(1)*x1(:,n1)-d1(2)*x1(:,1);
x(:,[2:n1])=x(:,[3:2:n2-1])-d1(1)*x1(:,[1:n1-1])-d1(2)*x1(:,[2:n1]);
x(:,[n1+1:n2])=x1(:,[1:n1]);
% d;
x(:,n2)=x(:,n2)-p1(1)*x(:,n1)-p1(2)*x(:,1);
x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])-p1(1)*x(:,[1:n1-1])-p1(2)*x(:,[2:n1]);
% p;
x(:,n1,:)=x(:,n1)-d2(1)*x(:,n2)-d2(2)*x(:,n1+1);
x(:,[1:n1-1])=x(:,[1:n1-1])-d2(1)*x(:,[n1+1:n2-1])-d2(2)*x(:,[n1+2:n2]);
% d;
x(:,n1+1)=x(:,n1+1)-p2(1)*x(:,n1-1)-p2(2)*x(:,n1);
x(:,n1+2)=x(:,n1+2)-p2(1)*x(:,n1)-p2(2)*x(:,1);
x(:,[n1+3:n2])=x(:,[n1+3:n2])-p2(1)*x(:,[1:n1-2])-p2(2)*x(:,[2:n1-1]);
% 归一
x(:,[1:n1])=d3*x(:,[1:n1]);
x(:,[n1+1:n2])=p3*x(:,[n1+1:n2]); clear x1;
n2=n2/2; % 原大小
n1=n2/2; % 一半大小
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 6.提升法反变换 else
n2=N/(2.^(T-1)); % 分解最小子块维数
n1=n2/2;
for time=1:T; % 行反变换
% 去归一
x([1:n1],:)=x([1:n1],:)/p3;
x([n1+1:n2],:)=x([n1+1:n2],:)/d3; % 反p;
x(n1,:)=x(n1,:)-p2(2)*x(n1+1,:)-p2(1)*x(n1+2,:);
x(n1-1,:)=x(n1-1,:)-p2(2)*x(n2,:)-p2(1)*x(n1+1,:);
x([1:n1-2],:)=x([1:n1-2],:)-p2(2)*x([n1+2:n2-1],:)-p2(1)*x([n1+3:n2],:);
% 反d;
x(n1+1,:)=x(n1+1,:)-d2(2)*x(n1,:)-d2(1)*x(1,:);
x([n1+2:n2],:)=x([n1+2:n2],:)-d2(2)*x([1:n1-1],:)-d2(1)*x([2:n1],:);
% 反p;
x1(1,:)=x(1,:)-p1(2)*x(n2,:)-p1(1)*x(n1+1,:);
x1([2:n1],:)=x([2:n1],:)-p1(2)*x([n1+1:n2-1],:)-p1(1)*x([n1+2:n2],:);
% 反d;
x(n2,:)=x(n2,:)-d1(2)*x1(n1,:)-d1(1)*x1(1,:);
x([2:2:n2-2],:)=x([n1+1:n2-1],:)-d1(2)*x1([1:n1-1],:)-d1(1)*x1([2:n1],:);
% 偶数
x([1:2:n2-1],:)=x1([1:n1],:);
clear x1;
% 列反变换
% 归一
x(:,[1:n1])=x(:,[1:n1])/d3;
x(:,[n1+1:n2])=x(:,[n1+1:n2])/p3; % 反d;
x(:,n1+1)=x(:,n1+1)+p2(1)*x(:,n1-1)+p2(2)*x(:,n1);
x(:,n1+2)=x(:,n1+2)+p2(1)*x(:,n1)+p2(2)*x(:,1);
x(:,[n1+3:n2])=x(:,[n1+3:n2])+p2(1)*x(:,[1:n1-2])+p2(2)*x(:,[2:n1-1]);
% 反p;
x(:,n1,:)=x(:,n1)+d2(1)*x(:,n2)+d2(2)*x(:,n1+1);
x(:,[1:n1-1])=x(:,[1:n1-1])+d2(1)*x(:,[n1+1:n2-1])+d2(2)*x(:,[n1+2:n2]);
% 反d;
x(:,n2)=x(:,n2)+p1(1)*x(:,n1)+p1(2)*x(:,1);
x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])+p1(1)*x(:,[1:n1-1])+p1(2)*x(:,[2:n1]);
% 反p;
x1(:,1)=x(:,1)+d1(1)*x(:,n2)+d1(2)*x(:,n1+1);
x1(:,[2:n1])=x(:,[2:n1])+d1(1)*x(:,[n1+1:n2-1])+d1(2)*x(:,[n1+2:n2]); % 奇偶
x(:,[2:2:n2])=x(:,[n1+1:n2]);
x(:,[1:2:n2-1])=x1(:,[1:n1]); clear x1;
n2=n2*2; % 原大小
n1=n2/2; % 一半大小 end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 7.结果输出 y=x;
% 传输最后结果 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 8.内存清理 clear x;
clear flag_max;
clear layer;
clear flag_trans;
clear N;
clear n1;
clear n2;
clear s;
clear ss;
clear u;
clear d1;
clear d2;
clear d3;
clear p1;
clear p2;
clear p3;
clear sa;
clear sb;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[ 本帖最后由 yejet 于 2006-8-31 20:32 编辑 ]
simon21
新科状元
精华: 6
积分: 400
帖子: 670
威望: 321 点
振动币: 23552 个
管理积分: 0 点
人气指数: 0
阅读权限: 80
注册: 2005-7-24
第 3 楼发表于 2006-3-31 08:38 资料 个人空间 短消息
2代小波示意程序
[Copy to clipboard]CODE:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 此程序用提升法实现第二代小波变换
%% 我用的是非整数阶小波变换
%% 采用时域实现,步骤先列后行
%% 正变换:分裂,预测,更新;
%% 反变换:更新,预测,合并
%% 只做一层(可以多层,而且每层的预测和更新方程不同) clear;clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
% 1.调原始图像矩阵 load wbarb; % 下载图像
f=X; % 原始图像
% f=[0 0 0 0 0 0 0 0 ;...
% 0 0 0 1 1 0 0 0 ;...
% 0 0 2 4 4 2 0 0 ;...
% 0 1 4 8 8 4 1 0 ;...
% 0 1 4 8 8 4 1 0 ;...
% 0 0 2 4 4 2 0 0 ;...
% 0 0 0 1 1 0 0 0 ;...
% 0 0 0 0 0 0 0 0 ;]; % 原始图像矩阵 N=length(f); % 图像维数
T=N/2;
% 子图像维数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%正变换%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 1.列变换
% A.分裂(奇偶分开) f1=f([1:2:N-1],:); % 奇数
f2=f([2:2:N],:); % 偶数 % f1(:,T+1)=f1(:,1); % 补列
% f2(T+1,:)=f2(1,:); % 补行 % B.预测 for i_hc=1:T;
high_frequency_column(i_hc,:)=f1(i_hc,:)-f2(i_hc,:);
end; % high_frequency_column(T+1,:)=high_frequency_column(1,:); % 补行 % C.更新 for i_lc=1:T;
low_frequency_column(i_lc,:)=f2(i_lc,:)+1/2*high_frequency_column(i_lc,:);
end; % D.合并
f_column([1:1:T],:)=low_frequency_column([1:T],:);
f_column([T+1:1:N],:)=high_frequency_column([1:T],:);
figure(1)
colormap(map);
image(f); figure(2)
colormap(map);
image(f_column);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 2.行变换
% A.分裂(奇偶分开) f1=f_column(:,[1:2:N-1]); % 奇数
f2=f_column(:,[2:2:N]); % 偶数
% f2(:,T+1)=f2(:,1); % 补行 % B.预测 for i_hr=1:T;
high_frequency_row(:,i_hr)=f1(:,i_hr)-f2(:,i_hr);
end; % high_frequency_row(:,T+1)=high_frequency_row(:,1); % 补行 % C.更新 for i_lr=1:T;
low_frequency_row(:,i_lr)=f2(:,i_lr)+1/2*high_frequency_row(:,i_lr);
end; % D.合并
f_row(:,[1:1:T])=low_frequency_row(:,[1:T]);
f_row(:,[T+1:1:N])=high_frequency_row(:,[1:T]);
figure(3)
colormap(map);
image(f_row);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%反变换%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 1.行变换
% A.提取(低频高频分开) f1=f_row(:,[T+1:1:N]); % 奇数
f2=f_row(:,[1:1:T]); % 偶数
% f2(:,T+1)=f2(:,1); % 补行 % B.更新 for i_lr=1:T;
low_frequency_row(:,i_lr)=f2(:,i_lr)-1/2*f1(:,i_lr);
end; % C.预测 for i_hr=1:T;
high_frequency_row(:,i_hr)=f1(:,i_hr)+low_frequency_row(:,i_hr);
end; % high_frequency_row(:,T+1)=high_frequency_row(:,1); % 补行
% D.合并(奇偶分开合并)
f_row(:,[2:2:N])=low_frequency_row(:,[1:T]);
f_row(:,[1:2:N-1])=high_frequency_row(:,[1:T]);
figure(4)
colormap(map);
image(f_row);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 2.列变换
% A.提取(低频高频分开) f1=f_row([T+1:1:N],:); % 奇数
f2=f_row([1:1:T],:); % 偶数 % f1(:,T+1)=f1(:,1); % 补列
% f2(T+1,:)=f2(1,:); % 补行 % B.更新 for i_lc=1:T;
low_frequency_column(i_lc,:)=f2(i_lc,:)-1/2*f1(i_lc,:);
end; % C.预测 for i_hc=1:T;
high_frequency_column(i_hc,:)=f1(i_hc,:)+low_frequency_column(i_hc,:);
end; % high_frequency_column(T+1,:)=high_frequency_column(1,:); % 补行 % D.合并(奇偶分开合并)
f_column([2:2:N],:)=low_frequency_column([1:T],:);
f_column([1:2:N-1],:)=high_frequency_column([1:T],:);
figure(5)
colormap(map);
image(f_column);
[ 本帖最后由 yejet 于 2006-8-31 20:32 编辑 ]
simon21
新科状元
精华: 6
积分: 400
帖子: 670
威望: 321 点
振动币: 23552 个
管理积分: 0 点
人气指数: 0
阅读权限: 80
注册: 2005-7-24
第 4 楼发表于 2006-3-31 08:40 资料 个人空间 短消息
回复:(simon21)[分享]个人收集的一些关于小波分析的...
二代小波漫谈
现在我就举例,对一个8点序列,怎样实现第二代小波变换。
1. 奇偶分开。
非常简单,就是[2,4,6,8]组成一列向量,[1,3,5,7]组成一列向量。
2. 预测。
用[2,4,6,8]来预测[1,3,5,7]。比如说1,3估计2; 3,5估计4; 5,7估计6; 7,1估计8。(边缘处理,我采用循环方法)。估计公式可以用别人的,也可以自己做。举一个线性的例子:2=1*a+3*b,4=3*a+5*b,...,其他的都一样。这样我们就可找到最优的a,b,使得(2-(1*a+3*b)).^2+(4-(3*a+5*b)).^2+...最小化。就是最小均方准则。若正好为零,说明偶可以完全预测奇,也就是我们只要存储偶数列向量,和a,b就可以了,压缩也就是实现了。对于信号很长序列,就等于压缩了一半。当然,我们可以采用更复杂的立方差值预测,多项式预测,或其它的准则,来使其最小,这样我们的压缩也就得到了最优。
3. 提升。
我们总希望,均方为零,但可望不可及。于是,提升就需要了。我们经过预测后,要存储的是偶数序列[2,4,6,8],新的奇数序列[n1,n3,n5,n7]=[2-(1*a+3*b),4-(3*a+5*b),...]和线性变换系数(a,b)。这里新的奇数序列就是高频分量。但偶数序列是不能完全代表信号的性质的,有所差距。所以我们要对偶数序列进行修正。即所谓的提升。我们这次用个简单的提升吧。[n2,n4,n6,n8]=[2,4,6,8]+k*[n1,n3,n5,n7]。[n2,n4,n6,n8],就是要分解的低频分量。那k怎么求呢?因为要保持n2,n4,n6,n8和原始信号[1,2,3,4,5,6,7,8]一样的性质。一般就是均值和高阶矩。这里就一个未知数k,所以用均值相等,就行了。1/8*(1+2+3+..8)=1/4*(n2+n4+n6+n8)。k很容易就求出来了。我们最终存储的就是[n1,n3,n5,n7]和[n2,n4,n6,n8]以及a,b,k。
现在,所谓的第二代就完了。再说几句。
1.反变换,就是3->2->1。
2.二维。先行提升,再列提升。(我置顶的贴子里有harr二维提升的源代码)。
3.整数阶。就是加一个取整。
4.多层或小波包提升,就是在对序列[n1,n3,n5,n7]或[n2,n4,n6,n8],再做1->2->3。
5.灵活。不一定是a,b,也可能就一个a,或a,b,c;不一定是一个k,也可能是k1,k2。但越多计算量太大。最好是用大师们做好的CDF,5/3,7/9等。
6.最重要的,任何一代小波,总可以通过一次或多次提升实现。它和一代小波没有本质区别。
7.优势。文献都有,我随便谈谈。时域实现,最优压缩,无边缘效应,灵活多变,无损压缩,编程方便,速度快。
文章写完了,希望对大家有帮助。最主要的,动手编,不要依赖MATLABM,这样才有所体会。希望和大家多交流。
simon21
新科状元
精华: 6
积分: 400
帖子: 670
威望: 321 点
振动币: 23552 个
管理积分: 0 点
人气指数: 0
阅读权限: 80
注册: 2005-7-24
第 5 楼发表于 2006-3-31 08:45 资料 个人空间 短消息
回复:(simon21)[分享]个人收集的一些关于小波分析的...
小波滤波器构造和消噪程序(2个)
1.重构
% mallet_wavelet.m
% 此函数用于研究Mallet算法及滤波器设计
% 此函数仅用于消噪
a=pi/8; %角度赋初值
b=pi/8;
%低通重构FIR滤波器h0(n)冲激响应赋值
h0=cos(a)*cos(b);
h1=sin(a)*cos(b);
h2=-sin(a)*sin(b);
h3=cos(a)*sin(b);
low_construct=[h0,h1,h2,h3];
L_fre=4; %滤波器长度
low_decompose=low_construct(end:-1:1); %确定h0(-n),低通分解滤波器
for i_high=1:L_fre; %确定h1(n)=(-1)^n,高通重建滤波器
if(mod(i_high,2)==0);
coefficient=-1;
else
coefficient=1;
end
high_construct(1,i_high)=low_decompose(1,i_high)*coefficient;
end
high_decompose=high_construct(end:-1:1); %高通分解滤波器h1(-n)
L_signal=100; %信号长度
n=1:L_signal; %信号赋值
f=10;
t=0.001;
y=10*cos(2*pi*50*n*t).*exp(-20*n*t);
figure(1);
plot(y);
title('原信号');
check1=sum(high_decompose); %h0(n)性质校验
check2=sum(low_decompose);
check3=norm(high_decompose);
check4=norm(low_decompose);
l_fre=conv(y,low_decompose); %卷积
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -