%%%%%%%%%%%%%%%%%%%%%%%%%%%ntvc xjh  2007-10-12%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%實驗二 上機內容(1)參考程序1：電子051、053班使用%%%%%%%%%%%%%%%%%%%%%%
%采用矩陣運算來計算頻率響應
w=[0:1:500]*pi/500;   %[0,pi]分為501個點
X1=1-exp(-1*j*w*6);
X2=1-exp(-1*j*w);
X22=X2+(X2==0)*eps;    %邏輯數組參加運算，使“0”被機器零代替
X=X1./X22;             %計算出頻率響應
magX=abs(X);            %取其幅度
angX=angle(X).*180./pi; %取其相位；
figure(1)
subplot(2,1,1);
plot(w/pi,magX);
title('幅頻響應');ylabel('幅度');
subplot(2,1,2);
plot(w/pi,angX);
title('相頻特性');ylabel('相位(以“度”為單位)');xlabel('以pi為單位');

%以下程序僅是為了調用zplane和impz函數來作零極點圖和單位脈沖響應
b=[1,0,0,0,0,0,-1];    %b是由分子多項式系數構成的數組；
a=[1,-1];              %a是由分母多項式系數構成的數組；
figure(2)
subplot(2,1,1);zplane(b,a);title('零極點圖');
subplot(2,1,2);impz(b,a);

%%%%%%%%%%%%%%%%%%%%%%%%上機內容(1)參考程序2：電子052班使用%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%使用MATLAB信號處理工具包中的函數freqz來計算頻率響應。b是由分子多項式系數構成的數組；
%a是由分母多項式系數構成的數組；
b=[1,0,0,0,0,0,-1];      %b是由分子多項式系數構成的數組；  
a=[1,-1];                %a是由分母多項式系數構成的數組；
[h,w]=freqz(b,a);      
am=20*log10(abs(h));     %%取其幅度,并轉換成以(dB)為單位
ph=angle(h);             %%取其相位。
w1=w/pi;                  %方便繪圖時橫軸以pi為單位
figure(1)
subplot(2,1,1);
plot(w1,am);title('幅頻響應');ylabel('幅度(dB)');
subplot(2,1,2);
plot(w1,ph);title('相頻特性');ylabel('相位(以“弧度”為單位)');xlabel('以pi為單位');
figure(2)
subplot(2,1,1);zplane(b,a);title('零極點圖');
subplot(2,1,2);impz(b,a);

%%%%%%%%%%%%%%%%%%%%%%%%%上機內容(2)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=[1,1,1,0,0,0,0];N=7;
nxtide=0:4*7-1;
xtide=x(mod(nxtide,N)+1);
stem(nxtide,xtide);axis([0 27 0 1]);

XK=dfs(x,N);                %dfs函數先自建，具體內容見附錄
magXK=abs(XK);
angXK=angle(XK).*180./pi;
K=[0:6];
subplot(2,1,1);stem(K,magXK);
title('幅度部分');ylabel('振幅');
subplot(2,1,2);stem(K,angXK);xlabel('K');ylabel('相位');
title('相位部分');axis([0 6 -180 180]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%附dfs函數內容,點擊MATLAB菜單File/new/M-file,將下列內容輸入后，并保存為函數名。

function  [XK]=dfs(x,N)
%計算離散傅里葉級數
%調用形式：[XK]=dfs(x,N)
%x是周期序列的一個周期序列值，0<=n<=N-1
%XK是DFS結果，0<=k<=N-1
%N是一個周期的長度
n=[0:1:N-1];
k=[0:1:N-1];
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
XK=x*WNnk;

%%%%%%%%%%%%%%%%%%%%%%%%%上機內容(3)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=[3,11,7,0,-1,4,2]; nx=[-3:3];
h=[2,3,0,-5,2,1];    nh=[-1:4];
[y,ny]=conv_m(x,nx,h,nh);             %先建conv_m函數
subplot(3,1,1);stem(nx,x);axis([-4 7 -1 11]);title('x(n)');
subplot(3,1,2);stem(nh,h);axis([-4 7 -5 3]);title('h(n)');
subplot(3,1,3);stem(ny,y);axis([-4 7 -51 50]);title('卷積結果y(n)=x(n)*h(n)');

%%%%%%%%%%%%%%%%%%%%%%%%%
%附conv_m函數內容

function [y,ny]=conv_m(x,nx,h,nh)
%conv進一步完善程序，可用來計算任意位置序列的卷積。
%調用格式：[y,ny]=conv_m(x,nx,h,nh)
%[y,ny]是卷積結果
%[x,nx]是第一個序列信號
%[h,nh]是第二個序列信號
%
nyb=nx(1)+nh(1);
nye=nx(length(x))+nh(length(h));
ny=[nyb:nye];
y=conv(x,h);


%%%%%%%%%%%%%%%%%%%%%%%%%上機內容(4)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(1,1,1);
n=-5:5;
x=(-0.9).^n;
k=-200:200;
w=(pi/100)*k;
X=x*(exp(-j*pi/100)).^(n'*k);      %求DTFT
magX=abs(X);                       %取幅度
angX=angle(X);                     %取相角

subplot(2,1,1);plot(w/pi,magX);grid
axis([-2 2 0 15]);
xlabel('頻率以pi為單位');ylabel('幅度|X|');

subplot(2,1,2);plot(w/pi,angX/pi);grid
axis([-2 2 -1 1]);
xlabel('頻率以pi為單位');ylabel('相位(rad/pi)');




%%%%%%%%%%%%%%%%%%%%%%%%%上機內容(5)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%附dft函數內容

function [Xk]=dft(x,N)
%該函數用來計算DFT
%調用形式：[XK]=dft(x,N)
%x是N點有限長序列值，0<=n<=N-1
%XK是DFT結果，0<=k<=N-1
%N是計算DFT的點數

n=[0:1:N-1];
k=[0:1:N-1];
WN=exp(-j*2*pi/N);   %旋轉因子
nk=n'*k;
WNnk=WN.^nk;
Xk=x*WNnk;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K1=[0:3];
x=[1 1 1 1];
N1=4;
Xk=dft(x,N1);
magXk=abs(Xk);              %DFT的模
phaXk=angle(Xk).*180./pi;   %DFT的相角

w=[0:1:500]*2*pi/500;
X1=1-exp(-j*w*4);
X2=1-exp(-j*w);
X22=X2+(X2==0)*eps;
X=X1./X22;
magX=abs(X);                %X連續傅里葉變換的模
K2=[0:11];
yn=[1 1 1 1 0 0 0 0 0 0 0 0 ];
N2=12;
Yk=dft(yn,N2);
magYk=abs(Yk);             %DFT的模
w=[0:1:500]*2*pi/500;
Y1=1+exp(-j*w)+exp(-j*w*2)+exp(-j*w*3);  %分子多項式
Y2=1;
Y=Y1./Y2;
magY=abs(Y);                %Y連續傅里葉變換的模

%以下是繪圖語句
subplot(2,2,1);stem(K1,x);title('信號x(n)');axis([0 4 -0.1 1.2]);
subplot(2,2,2);plot(2*w/pi,magX,'--');hold on;title('X幅度函數及抽樣');
ylabel('振幅');stem(K1,magXk);axis([0 4 -0.5 5]);
subplot(2,2,3); stem(K2,yn);title('信號y(n)');axis([0,12,-0.1,1.2]);
subplot(2,2,4); plot(6*w/pi,magY,'--');hold on;
title('Y幅度函數及抽樣');ylabel('振幅');
stem(K2,magYk);axis([0,12,-0.5,5]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%