样条小波的双尺度关系

题目：验证样条小波的双尺度关系，绘制样条小波的尺度函数和小波函数%%构建尺度函数functionmea_w=funmea_w(m,w,N)mea_w=(2*sin(w/2)./(w+eps)).^m.*exp(-j*m*w/2);%求Nm(w)%求sum（|N2m(w+2kpi)|^2）,利用循环mea_w1=N(2*m,m+1);fori=-m+1:-1mea_w1=mea_w1+2*N(2*m,m+i+1)*cos(i*w);endmea_w=mea_w./sqrt(mea_w1);%求出尺度函数end%%构建插值函数functionN=nm(m,k)N=zeros(m,k);N(1,1)=1;ifm1fori=2:mforj=2:kN(i,j)=(j-1)*N(i-1,j)/(i-1)+(i-j+1)*N(i-1,j-1)/(i-1);%利用插值公式计算n2一直到nm。0时刻对应的值为数组素引1，其他依次类推endendendend%%主程序m=4;N=nm(2*m,2*m+1);%求样条函数在整数点的值,利用查值公式w=0:0.1*pi:10*pi;mea_w=funmea_w(m,w,N);%求尺度函数figure(1);subplot(3,2,1),plot(w/pi,abs(mea_w));%画出尺度函数的幅频特性title('尺度函数的幅频图(4阶)');mea_ang=atan(imag(mea_w)./(real(mea_w)+eps));%计算相位角subplot(3,2,3),plot(w/pi,mea_ang);%画出尺度函数的相频特性title('尺度函数的相频图(4阶)');mea_x=ifft(mea_w);%求反傅立叶变换，求出尺度函数的时域函数。mea_x=ifftshift(mea_x);%对时域函数进行做一些搬移，为了放在中间。方便观察subplot(3,2,5),plot(w/pi,real(mea_x));%画出尺度函数的时域图title('尺度函数的时域图(4阶)');mea_2w=funmea_w(m,2*w,N);Hw=mea_2w./(mea_w+eps);%求H(w),利用H_k=ifft(Hw(1:21));%求h（k）利用H_kH_kreal=real(H_k)wdiv2=w/2;Hwpidiv2=funmea_w(m,2*(wdiv2+pi),N)./funmea_w(m,wdiv2+pi,N);%求H(w+pi)Hwpidiv2_=real(Hwpidiv2)-i*imag(Hwpidiv2);%求H(w/2+pi)的共轭G_wdiv2=exp(-j*wdiv2).*Hwpidiv2_;%求G(w/2)mea_wdiv2=funmea_w(m,wdiv2,N);%求尺度函数（w/2）bwav_w=G_wdiv2.*mea_wdiv2;%求小波函数figure(2);subplot(3,2,1),plot(w/pi,abs(bwav_w));%画出小波函数的幅频特性title('小波函数的幅频图(4阶)');bwav_ang=atan(imag(bwav_w)./(real(bwav_w)+eps));%计算相位角subplot(3,2,3),plot(w/pi,bwav_ang);%画出小波函数的相频特性title('小波函数的相频图(4阶)');bwav_x=ifft(bwav_w);%求反傅立叶变换，求出小波函数的时域函数。bwav_x=ifftshift(bwav_x);%对时域函数进行做一些搬移，为了放在中间。方便观察subplot(3,2,5),plot(w/pi,real(bwav_x));%画出小波函数的时域图title('小波函数的时域图(4阶)');%%m=5m=5;N=nm(2*m,2*m+1);%求样条函数在整数点的值,利用查值公式w=0:0.1*pi:10*pi;mea_w=funmea_w(m,w,N);%求尺度函数figure(1);subplot(3,2,2),plot(w/pi,abs(mea_w));%画出尺度函数的幅频特性title('尺度函数的幅频图(5阶)');mea_ang=atan(imag(mea_w)./(real(mea_w)+eps));%计算相位角subplot(3,2,4),plot(w/pi,mea_ang);%画出尺度函数的相频特性title('尺度函数的相频图(5阶)');mea_x=ifft(mea_w);%求反傅立叶变换，求出尺度函数的时域函数。mea_x=ifftshift(mea_x);%对时域函数进行做一些搬移，为了放在中间。方便观察subplot(3,2,6),plot(w/pi,real(mea_x));%画出尺度函数的时域图title('尺度函数的时域图(5阶)');mea_2w=funmea_w(m,2*w,N);Hw=mea_2w./(mea_w+eps);%求H(w),利用H_k=ifft(Hw(1:21));%求h（k）利用H_kH_kreal=real(H_k)wdiv2=w/2;Hwpidiv2=funmea_w(m,2*(wdiv2+pi),N)./funmea_w(m,wdiv2+pi,N);%求H(w+pi)Hwpidiv2_=real(Hwpidiv2)-i*imag(Hwpidiv2);%求H(w/2+pi)的共轭G_wdiv2=exp(-j*wdiv2).*Hwpidiv2_;%求G(w/2)mea_wdiv2=funmea_w(m,wdiv2,N);%求尺度函数（w/2）bwav_w=G_wdiv2.*mea_wdiv2;%求小波函数figure(2);subplot(3,2,2),plot(w/pi,abs(bwav_w));%画出小波函数的幅频特性title('小波函数的幅频图(5阶)');bwav_ang=atan(imag(bwav_w)./(real(bwav_w)+eps));%计算相位角subplot(3,2,4),plot(w/pi,bwav_ang);%画出小波函数的相频特性title('小波函数的相频图(5阶)');bwav_x=ifft(bwav_w);%求反傅立叶变换，求出小波函数的时域函数。bwav_x=ifftshift(bwav_x);%对时域函数进行做一些搬移，为了放在中间。方便观察subplot(3,2,6),plot(w/pi,real(bwav_x));%画出小波函数的时域图title('小波函数的时域图(5阶)');051000.51尺度函数的幅频图(4阶)0510-202尺度函数的相频图(4阶)0510-0.200.2尺度函数的时域图(4阶)051000.51尺度函数的幅频图(5阶)0510-202尺度函数的相频图(5阶)0510-0.100.1尺度函数的时域图(5阶)051000.51小波函数的幅频图(4阶)0510-202小波函数的相频图(4阶)0510-0.200.2小波函数的时域图(4阶)051000.51小波函数的幅频图(5阶)0510-202小波函数的相频图(5阶)0510-0.200.2小波函数的时域图(5阶)