wavelet_char.m

小波压缩编码平台软件包

function [f1freq,f1phase,f2freq,f2phase]=wavelet_char(wavelet1,wavelet2,type,plt);
%[f1freq,f1phase,f2freq,f2phase] = wavelet_char(wavelet1,wavelet2,type,plt);
%Version: 3.10, Date: 2006/01/31, author: Nikola Sprljan
%Computes (and plots) frequency and phase charateristic of a wavelet
%
%Input: 
% wavelet1 - identification string of the first wavelet 
% wavelet2 - identification string of the second wavelet 
% type - specifies wavelet type
%         if ('d' in type) analysis (decomposition) filter
%         if ('r' in type) synthesis (reconstruction) filter
%         if ('l' in type) then low-pass
%         if ('h' in type) then high-pass
% plt - [optional, default = ''] 
%        if (plt == 'plot') then plots the characteristics
%
%Output: 
% f1freq,f2freq - frequency charateristics
% f1phase,f2phase - phase charateristics
%
%Uses: 
% load_wavelet.m (Wavelet Toolbox)
%
%Note: 
% Number of discrete frequencies in which transfer function is computed is 
% specified by wsamples variable
%
%Example:
% wavelet_char('CDF_9x7','Haar','dl','plot');
% [f1f,f1p]=wavelet_char('haar','haar','d');

wsamples=1024; %change this to change sampling of the frequency axis
if nargin<=3
    plt='';
end;
if strcmp(wavelet1,wavelet2) %if it is the same wavelet
    wvf=load_wavelet(wavelet1);
    switch type(1)
        case 'd'
            filter1 = wvf.filt_AL;
            leg1=' low-pass';
            filter2 = wvf.filt_AH;
            leg2=' high-pass';
        case 'r'
            filter1 = wvf.filt_SL;
            leg1=' low-pass';
            filter2 = wvf.filt_SH;
            leg2=' high-pass';
        case 'l'
            filter1 = wvf.filt_AL;
            leg1=' analysis';
            filter2 = wvf.filt_SL;
            leg2=' synthesis';
        case 'h'
            filter1 = wvf.filt_AH;
            leg1=' analysis';
            filter2 = wvf.filt_SH;
            leg2=' synthesis';
    end;
else
    wvf1=load_wavelet(wavelet1);
    wvf2=load_wavelet(wavelet2);
    switch type
        case {'dl','ld'}
            filter1 = wvf1.filt_AL;
            leg1=' analysis low-pass';
            filter2 = wvf2.filt_AL;
            leg2=' analysis low-pass';
        case {'dh','hd'}
            filter1 = wvf1.filt_AH;
            leg1=' analysis high-pass';
            filter2 = wvf2.filt_AH;
            leg2=' analysis high-pass';
        case {'rl','lr'}
            filter1 = wvf1.filt_SL;
            leg1=' sytnhesis low-pass';
            filter2 = wvf2.filt_SL;
            leg2=' sytnhesis low-pass';
        case {'rh','hr'}
            filter1 = wvf1.filt_SH;
            leg1=' sytnhesis high-pass';
            filter2 = wvf2.filt_SH;
            leg2=' analysis high-pass';
    end;
end;

%%computed by FFT 
%filter1
% f1freq = abs(fft(filter1,wsamples*2));
% f1freq = f1freq(1:size(f1freq,2)/2);
% %filter2
% f2freq = abs(fft(filter2,wsamples*2));
% f2freq = f2freq(1:size(f2freq,2)/2);
% % h = fir1(20,1/2); %FIR example leftover
% % W3=abs(fft(h,wsamples*2)); %FIR example leftover
% % W3 = W3(1:size(W3,2)/2);

w = 0:pi/(wsamples-1):pi;
%%computed by z-transform
%filter1
W = repmat(w,[length(filter1) 1]);
z = cos(W)+i*sin(W);
n = 0:-1:-length(filter1)+1;
N = repmat(n',[1 wsamples]);
F1 = repmat(filter1',[1 wsamples]);
F1 = sum(F1.*(z.^N));
f1freq=abs(F1);
f1phase=angle(F1);
%filter2
W = repmat(w,[length(filter2) 1]);
z = cos(W)+i*sin(W);
n = 0:-1:-length(filter2)+1;
N = repmat(n',[1 wsamples]);
F2 = repmat(filter2',[1 wsamples]);
F2 = sum(F2.*(z.^N)); 
f2freq=abs(F2);
f2phase=angle(F2); 

if strcmp(plt,'plot')
 wavelet1=strrep(wavelet1,'_','');
 wavelet2=strrep(wavelet2,'_',''); 
 %%%frequency plot%%%
 maxf1freq=max(f1freq(:));
 maxf2freq=max(f2freq(:));
 maxyos=1.1*max([maxf1freq maxf2freq]);

 %figure('Name','Wavelet filter fft spectrum'); %FFT
 figure('Name',['Wavelet filter z-transfer function (' wavelet1 ',' wavelet2 ')'],'NumberTitle','off'); %z-transform
 plot(w,f1freq,'-k',w,f2freq,'-r','LineWidth',2);
 %plot(x,f1freq,'-y',x,f2freq,'-k',x,W3(,'-b','LineWidth',2); %FIR example leftover
 
 %y axis
 text(-0.4,maxyos/2,'|A(\omega)|','FontSize',20,'Rotation',90,'HorizontalAlignment','center');
 %x axis
 text(pi/2,-maxyos/30,'  0                 \pi/4                \pi/2                 3\pi/4                  \pi',...
     'FontSize',14,'HorizontalAlignment','center');
 text(pi/2,-maxyos/10,'\omega','FontSize',20,'HorizontalAlignment','center');
 
 axis([0 pi 0 maxyos]); 
 %set(gcf,'Color',[1 1 1]);
 set(gca,'FontSize',15,'XTick',0:pi/4:pi,'LineWidth',2);
 set(gca,'XTickLabel',{});
 
 legend([wavelet1 leg1],[wavelet2 leg2],0);
 grid on;
 %set(gca,'Color','none');
 %set(gcf,'Color','none');
 
 %%%phase plot%%%
 figure('Name',['Wavelet filter phase characteristic (' wavelet1 ',' wavelet2 ')'],'NumberTitle','off');
 plot(w,f1phase,'-k',w,f2phase,'-r','LineWidth',2);
 %y axis
 prop_name(1)={'FontSize'};
 properti(1)={14};
 prop_name(2)={'HorizontalAlignment'};
 properti(2)={'right'};
 text(-0.1,pi,'\pi',prop_name,properti);
 text(-0.1,pi/2,'\pi/2',prop_name,properti);
 text(-0.1,0,'0',prop_name,properti);
 text(-0.1,-pi/2,'-\pi/2',prop_name,properti);
 text(-0.1,-pi,'-\pi',prop_name,properti);
 text(-0.4,0,'\Phi(\omega)','FontSize',20,'Rotation',90,'HorizontalAlignment','center');
 %x axis 
 text(pi/2,-pi-2*pi/30,'  0                 \pi/4                \pi/2                 3\pi/4                  \pi',...
     'FontSize',14,'HorizontalAlignment','center');
 text(pi/2,-pi-2*pi/10,'\omega','FontSize',20,'HorizontalAlignment','center');
 
 axis([0 pi -pi pi]); 
 %set(gcf,'Color',[1 1 1]);
 set(gca,'FontSize',15,'XTick',0:pi/4:pi,'YTick',-pi:pi/2:pi,'LineWidth',2);
 set(gca,'XTickLabel',{});
 set(gca,'YTickLabel',{});
 
 legend([wavelet1 leg1],[wavelet2 leg2],0);
 grid on;
 %set(gca,'Color','none');
 %set(gcf,'Color','none');
end;