fdkaiser.m

matlabDigitalSigalProcess内有文件若干

            end
        end
    else
        xd = L3_2.xdata;
        newFc2 = inbounds(xd(1),[Wn2.range(1)+minspacing Fs/2-minspacing]);
        if newFc2 ~= Wn2.value
            Wn2.value = newFc2;
            Wn1.range = [0 newFc2];
            if newFc2 ~= xd(1)
                L3_2.xdata = newFc2([1 1]);
            end
        end
    end

%---------------------------------------------------------------------
% fdkaiser('L1drag',type,ind)
%  vertex drag callback of L1 - passband line
%  Inputs:
%     type - band configuration 1==low, 2=high, 3=pass, 4=stop
%     ind - index of vertex being dragged
case 'L1drag'
    type = varargin{2};
    ind = varargin{3};
    xd = L1.xdata;
    minspacing = Fs/500;
    switch type
    case 1  % lowpass
        newf1 = inbounds(xd(ind),[minspacing f2.value-minspacing]);
        xd([2 5]) = newf1;
        L1.xdata = xd;
        f1.value = newf1;
        i = 1;
    case 2  % highpass
        newf2 = inbounds(xd(ind),[f1.value+minspacing Fs/2-minspacing]);
        xd([1 4]) = newf2;
        L1.xdata = xd;
        f2.value = newf2;
        i = 2;
    case 3  % bandpass
    % L1 xdata = [f2 f3 NaN f2 f3]
        if any(ind == [1 4])   % dragging f2
            newf2 = inbounds(xd(ind),[f1.value+minspacing f3.value-minspacing]);
            xd([1 4]) = newf2;
            L1.xdata = xd;
            f2.value = newf2;
            i = 2;
        else % dragging f3
            newf3 = inbounds(xd(ind),[f2.value+minspacing f4.value-minspacing]);
            xd([2 5]) = newf3;
            L1.xdata = xd;
            f3.value = newf3;
            i = 3;
        end
    case 4   % bandstop
    % L1 xdata = [0 f1 NaN 0 f1 NaN f4 Fs/2 NaN f4 Fs/2]
        if any(ind == [2 5])   % dragging f1
            newf1 = inbounds(xd(ind),[minspacing f2.value-minspacing]);
            xd([2 5]) = newf1;
            L1.xdata = xd;
            f1.value = newf1;
            i = 1;
        else % dragging f4
            newf4 = inbounds(xd(ind),[f3.value+minspacing Fs/2-minspacing]);
            xd([7 10]) = newf4;
            L1.xdata = xd;
            f4.value = newf4;
            i = 4;
        end
    end
    if ~minOrdCheckbox.value
        fChangeMeas(i,type,Fs,Lresp,L1,L2,f1,f2,f3,f4,Rp,Rs)
        return
    end
    ax.xlimPassband = fdutil('xlimpassband',type,...
                         Fs,f1.value,f2.value,f3.value,f4.value);

%---------------------------------------------------------------------
% fdkaiser('L2drag',type,ind)
%  drag callback of L2 - stopband line
%  Inputs:
%     type - band configuration 1==low, 2=high, 3=pass, 4=stop
%     ind - index of vertex being dragged
case 'L2drag'
    type = varargin{2};
    ind = varargin{3};
    xd = L2.xdata;
    minspacing = Fs/500;
    switch type
    case 1  % lowpass
        newf2 = inbounds(xd(ind),[f1.value+minspacing Fs/2-minspacing]);
        xd(1) = newf2;
        L2.xdata = xd;
        f2.value = newf2;
        i = 2;
    case 2  % highpass
        newf1 = inbounds(xd(ind),[minspacing f2.value-minspacing]);
        xd(2) = newf1;
        L2.xdata = xd;
        f1.value = newf1;
        i = 1;
    case 3  % bandpass
    % L2 xdata = [0 f1 NaN f4 Fs/2]
        if ind == 2   % dragging f1
            newf1 = inbounds(xd(ind),[minspacing f2.value-minspacing]);
            xd(2) = newf1;
            L2.xdata = xd;
            f1.value = newf1;
            i = 1;
        else % dragging f4
            newf4 = inbounds(xd(ind),[f3.value+minspacing Fs/2-minspacing]);
            xd(4) = newf4;
            L2.xdata = xd;
            f4.value = newf4;
            i = 4;
        end
    case 4   % bandstop
    % L2 xdata = [f2 f3]
        if ind == 1   % dragging f2
            newf2 = inbounds(xd(ind),[f1.value+minspacing f3.value-minspacing]);
            xd(1) = newf2;
            L2.xdata = xd;
            f2.value = newf2;
            i = 2;
        else % dragging f3
            newf3 = inbounds(xd(ind),[f2.value+minspacing f4.value-minspacing]);
            xd(2) = newf3;
            L2.xdata = xd;
            f3.value = newf3;
            i = 3;
        end
    end
    if ~minOrdCheckbox.value
        fChangeMeas(i,type,Fs,Lresp,L1,L2,f1,f2,f3,f4,Rp,Rs)
    end
    
%---------------------------------------------------------------------
% fdkaiser('Rpdrag',type,ind)
%  drag callback of L1 - passband line
%  Inputs:
%     type - band configuration 1==low, 2=high, 3=pass, 4=stop
%     ind - index of segment being dragged
case 'Rpdrag'
    type = varargin{2};
    ind = varargin{3};
    yd = L1.ydata;
    
    switch ind
    case {4,10}    % dragging lower line
        below = yd(ind);
        if get(minOrdCheckbox,'value')==1    % estimate order
        % drag indicates change in desired spec for Rp 
            if below >= 0 
                below = -.00001;
            elseif below < -10
                below = -10;
            end
            dev = 1-10^(below/20);
            above = 20*log10(1+dev);
        else
        % drag indicates changing measurement
            above = yd(1);    
        end
        RpInd = 1;
    case {1,7}    % dragging upper line
        above = yd(ind);
        if get(minOrdCheckbox,'value')==1    % estimate order
        % drag indicates change in desired spec for Rp 
            if above > 10 
                above = 10;
            elseif above <= 0
                above = .00001;
            end
            dev = 10^(above/20)-1;
            below = 20*log10(1-dev);
        else
        % drag indicates changing measurement
            below = yd(4);    
        end
        RpInd = 2;
    end
    
    newRp = above-below;

    if minOrdCheckbox.value
        if type ~= 4
        % 'ydata':[maxpass maxpass NaN minpass minpass]
            yd = [above above NaN below below];
        else
        % 'ydata': [ maxpass maxpass NaN minpass minpass NaN ...
        %                    maxpass maxpass NaN minpass minpass])
            yd = [above above NaN below below ...
                   NaN above above NaN below below ];
        end
        L1.ydata = yd;
        Rp.value = newRp;
        ylim = [below above];
        dyl = (ylim(2)-ylim(1))*.15;
        ax.ylimPassband = ylim + [-dyl/2 dyl/2];
    else
        % Rp.value = newRp;
        RChangeMeas(1,type,Fs,Lresp,L1,L2,f1,f2,f3,f4,Rp,Rs,[below above],RpInd)
    end
 
%---------------------------------------------------------------------
% fdkaiser('Rsdrag',type,ind)
%  drag callback of L2 - stopband line
%  Inputs:
%     type - band configuration 1==low, 2=high, 3=pass, 4=stop
%     ind - index of segment being dragged
case 'Rsdrag'
    type = varargin{2};
    ind = varargin{3};
  
    yd = L2.ydata;
    newRs = -yd(ind);
    if newRs < 0
        newRs = 0;
    end
    if minOrdCheckbox.value
        switch type
        case {1,2,4}
            L2.ydata = [-newRs -newRs];
        case 3
            L2.ydata = [-newRs -newRs NaN -newRs -newRs];
        end
        set(Rs,'value', newRs)
    else
        set(Rs,'value', newRs)
        RChangeMeas(0,type,Fs,Lresp,L1,L2,f1,f2,f3,f4,Rp,Rs)
    end
    
end  % of function switch-yard

%---------------------------------------------------------------------
% -------- LOCAL FUNCTIONS START HERE  --------- 
%---------------------------------------------------------------------

function sethelp
%  following are common to several modules
global minOrdCheckbox bandpop order pbspecs sbspecs pbmeasures sbmeasures
global passframe stopframe passframe1 stopframe1
%  following static to fdkaiser
global f1 f2 f3 f4 Rp Rs Wn1 Wn2 Beta 
global Wn1m Wn2m Betam Rpm Rsm 
global order0 order1
global ax L1 L2 Lresp L3_1 L3_2
global Fs

% disp('setting help ... stub')


function [n,Wn,beta_val] = estimateOrder(type,Rp,Rs,Fs,f1,f2,f3,f4)
% [n,Wn,beta_val] = estimateOrder(type,Rp,Rs,Fs,f1,f2,f3,f4)
%   estimate filter order
%   takes the specifications as given by the input
%   parameters and estimates the order, and beta value
%   needed to meet those specifications.
%   Also returns Wn - cutoff frequencies for fir1 filter.
% Inputs:
%  type - 1,2,3,4 specifies band configuration
%  Rp, Rs passband, stopband ripple
%  Fs - sampling frequency
%  f1,f2  first two frequencies in ascending order
%  f3,f4 only needed if type == 3 or 4, remaining frequencies
%    f1,f2,f3,f4 are assumed between 0 and Fs/2 on input
% Outputs:
%  n - filter order
%  Wn - filter band edges for fir1, normalized to range [0...1]
%  beta_val - beta parameter of Kaiser window

% compute deviations and estimate order
if type == 1    % low pass
    dev = [ (10^(Rp/20)-1)/(10^(Rp/20)+1)  10^(-Rs/20) ];
    f = [f1 f2];  m = [1 0];
elseif type == 2   % high pass
    dev = [ 10^(-Rs/20) (10^(Rp/20)-1)/(10^(Rp/20)+1)];
    f = [f1 f2];  m = [0 1];
elseif type == 3   % band pass
    dev = [ 10^(-Rs/20) (10^(Rp/20)-1)/(10^(Rp/20)+1) 10^(-Rs/20)];
    f = [f1 f2 f3 f4];  m = [0 1 0];
elseif type == 4   % band stop
    dev = [(10^(Rp/20)-1)/(10^(Rp/20)+1) ...
              10^(-Rs/20) (10^(Rp/20)-1)/(10^(Rp/20)+1)];
    f = [f1 f2 f3 f4];  m = [1 0 1];
end
[n,Wn,beta_val] = kaiserord(f,m,dev,Fs);
Wn = Wn(:)';   % make Wn a row 
 

function yd = passbandlimits(Rp)
% return ydata = [minpass maxpass] of passband 
%   given Rp decibels of ripple
%   in passband (centered at 1 in linear scale)
    dev = (10^(Rp/20)-1)/(10^(Rp/20)+1);
    above = 20*log10(1+dev); below = 20*log10(1-dev);
    yd = [below above];

function [maxpass,minpass,minstop] = getMagMeasurements(ff,Hlog,type,...
                                          f1,f2,f3,f4);
%getMagMeasurements
% Finds passband and stopband ripple for given band edges 
% given a filter's magnitude response
% Inputs:
%  ff - xdata of response
%  Hlog - magnitude of response at frequencies ff, in dB
%  type - band configuration (lp = 1, hp = 2, bp = 3, bs = 4)
%  f1, f2, f3, f4 - band edges (f3 and f4 ignored if type < 3)
%           in same units as ff
% Output:
%  maxpass - maximum passband value - 2 element vector if type == 4 (bandstop)
%  minpass - minimum passband value - 2 element vector if type == 4 (bandstop)
%  minstop - maximum stopband value - 2 element vector is type == 3 (bandpass)
    switch type
    case 1   % lowpass
        passInd = find(ff<=f1);
        stopInd = find(ff>=f2);
    case 2   % highpass
        stopInd = find(ff<=f1);
        passInd = find(ff>=f2);
    case 3   % bandpass
        stopInd1 = find(ff<=f1);
        stopInd2 = find(ff>=f4);
        [Hmax,centerInd] = max(Hlog);
        passInd1 = find((ff>=f2)&(ff<=ff(centerInd)));
        passInd2 = find((ff>ff(centerInd))&(ff<=f3));
    case 4   % bandstop
        passInd1 = find(ff<=f1);
        passInd2 = find(ff>=f4);       
        [peakInd,peaks] = fdutil('findpeaks',Hlog);
        [Hmin,ind] = min(peaks);
        centerInd = peakInd(ind);        
        stopInd1 = find((ff>=f2)&(ff<=ff(centerInd)));
        stopInd2 = find((ff>ff(centerInd))&(ff<=f3));
    end
    switch type
    case {3,4}
        maxpass = [max(Hlog(passInd1)) max(Hlog(passInd2))];
        minpass = [min(Hlog(passInd1)) min(Hlog(passInd2))];
        minstop = [max(Hlog(stopInd1)) max(Hlog(stopInd2))];
    otherwise
        maxpass = max(Hlog(passInd));
        minpass = min(Hlog(passInd));
        minstop = max(Hlog(stopInd));
    end

    
function [fm,Rp_range,Rs] = getFreqMeasurements(passbandChange,ff,Hlog,type,Rp_range,...
                               Rs,RpInd)
%getFreqMeasurements
% Finds band edges for a given passband and stopband ripple
% given a filter's magnitude response
% Inputs:
%  passbandChange - 1 if given Rp is a desired quantity, 0 if Rp is known valid
%  ff - xdata of response (assumed a column vector)
%  Hlog - magnitude of response at frequencies ff, in dB
%               (assumed a column vector)
%  type - band configuration (lp = 1, hp = 2, bp = 3, bs = 4)
%  Rp_range - this is either a scalar Rp = passband ripple, in dB
%        or a 2 element range Rp_range = [below above], diff(Rp_range) = Rp 
%  Rs - stopband attenuation in dB
%  RpInd - index into Rp_range indicating which passband line is being dragged
%      1 ==> lower, 2 ==> upper, [] ==> neither.  Defaults to [].
% Output:
%  fm - 2 or 4 element frequency vector in ascending order
%  Rp, Rs - actual value of ripple across entire band
%  Rp_range - two element vector [below above], diff(Rp_range) = Rp 
    Hlog = Hlog(:);
    ff = ff(:);
    
    if nargin < 7
        RpInd = [];
    end
    stopInd = find(Hlog<=-Rs);
    if length(Rp_range)==1
        if passbandChange
            Rp = Rp_range;
            switch type 
            case 1