fgp852.m

digital signal processing-MATLAB

%《数字信号处理教程——MATLAB释义与实现》第八章图8.5.2生成程序fgp852
% 二阶全通数字频带变换函数的四组曲线
% 电子工业出版社出版  陈怀琛编著 2004年9月
%
clear,close all
w=-linspace(-pi,pi,501);                    % 自变量相角数组的设定
alpha=input('alpha=  ');                    % 设定或输入参数alpha1
if isempty(alpha) 
    alpha=0.3;
end
alpha2=input('alpha2=  ');                  % 设定或输入参数alpha2
if isempty(alpha2) 
    alpha2=0.5;
end
alpha1=alpha;
q=exp(-j*w);                                % z平面自变量q的设定
Q1=(q-alpha)./(1-alpha*q);                  % 用表8.5.2中第一行的变换公式求Z平面映射
angQ1=angle(Q1);                            % 求Z平面的映射的相角
set(gcf,'color','w')                        % 置图形背景色为白
figure(1),subplot(2,2,1)
np1=find(abs(angQ1)<=1);                    % Z平面通带（设为1）的下标
wp1=w(np1);angQ1p=angQ1(np1);               % 找到对应的z平面自变量范围wp
plot(-angQ1,w),hold on,grid on              % 画映射曲线
plot(-angQ1p,wp1,'.'),                      % 画通带映射曲线
plot(ones(size(wp1))*4,wp1,'.r'),           % 在右边界(x=4)画z平面通带范围
plot(-angQ1p,-4*ones(size(wp1)),'.g'),      % 在下边界(y=-4)画Z平面通带范围
ylabel('待求低通的数字频率')
title('(a)数字低通原型变换为数字低通')
line([-1,-1],[-pi,pi]),line([1,1],[pi,-pi]) % 画z平面通带投影线
subplot(2,2,2)
Q2=-(q+alpha)./(1+alpha*q);                 % 表8.5.2中第二行的变换公式求Z平面映射
angQ2=angle(Q2);                            % 求Z平面的映射的相角
np2=find(abs(angQ2)<=1);                    % Z平面通带（设为1）的下标
wp2=w(np2);angQ2p=angQ2(np2);               % 找到对应的z平面自变量范围wp
plot(-angQ2,w),hold on,grid on             % 画映射曲线
plot(-angQ2p,wp2,'.'),                      % 画通带映射曲线
plot(ones(size(wp2))*4,wp2,'.r'),           % 在右边界(x=4)画z平面通带范围
plot(-angQ2p,-4*ones(size(wp2)),'.g'),
ylabel('高通滤波器数字频率')
title('(b)数字低通原型变换为数字高通')
line([-1,-1],[-pi,pi]),line([1,1],[pi,-pi]) % 画z平面通带投影线
subplot(2,2,3)
Q3=-(q.^2-alpha1*q+alpha2)./(alpha2*q.^2-alpha1*q+1); % 表8.5.2中第三行的变换公式求Z平面映射
angQ3=unwrap(angle(Q3));                              % 求Z平面的映射的相角
np3=find(abs(angle(Q3))<=1);                           % Z平面通带（设为1）的下标
wp3=w(np3);angQ3p=angle(Q3(np3));                % 找到对应的z平面自变量范围wp
plot(-angle(Q3),w),hold on,grid on             % 画映射曲线
plot(-angQ3p,wp3,'.'),                      % 画通带映射曲线
plot(ones(size(wp3))*4,wp3,'.r'),           % 在右边界(x=4)画z平面通带范围
plot(-angQ3p,-4*ones(size(wp3)),'.g'),      % 在下边界(y=-4)画Z平面通带范围
xlabel('低通数字滤波器归一化频率(-pi~pi)'),
ylabel('带通滤波器数字频率')
title('(c)数字低通原型变换为数字带通')
line([-1,-1],[-pi,pi]),line([1,1],[pi,-pi]) % 画z平面通带投影线
subplot(2,2,4)
Q4=(q.^2-alpha1*q+alpha2)./(alpha2*q.^2-alpha1*q+1);  % 表8.5.2中第四行的变换公式求Z平面映射
angQ4=unwrap(angle(Q4)-2*pi);                          % 求Z平面的映射的相角
np4=find(abs(angle(Q4))<=1);                           % Z平面通带（设为1）的下标
wp4=w(np4);angQ4p=angle(Q4(np4));                        % 找到对应的z平面自变量范围wp
plot(-angle(Q4),w),hold on,grid on                   % 画映射曲线
plot(-angQ4p,wp4,'.'),                      % 在右边界(x=4)画z平面通带范围
plot(ones(size(wp4))*4,wp4,'.r'),           % 在下边界(y=-4)画Z平面通带范围
plot(-angQ4p,-4*ones(size(wp4)),'.g'),
xlabel('低通数字滤波器归一化频率(-pi~pi)'),
ylabel('带阻滤波器数字频率')
title('(d)数字低通原型变换为数字带阻')
line([-1,-1],[-pi,pi]),line([1,1],[pi,-pi]) % 画z平面通带投影线