% exa060202_cohere.m,  for example 6.2.2 
% to test cohere.m  ;

clear all;
% 構造低通濾波系數b1和高通濾波序列b2
Fs=1;N=1024;
A=[1 1  0 0];
f=[0 .6 .7 1];
weigh=[1 10];	
b1=remez(42,f,A,weigh);

A=[0 0  1 1];
f=[0 .5 .7 1];
weigh=[10 1];	
b2=remez(42,f,A,weigh);

[h1,w]=freqz(b1,1,256,1);
h1=abs(h1);
h1=20*log10(h1);
subplot(321);
plot(w,h1);
grid;
ylabel('低通濾波器')
axis([0 0.5 -70 10])

[h2,w]=freqz(b2,1,256,1);
h2=abs(h2);
h2=20*log10(h2);
subplot(322);
plot(w,h2);
grid;
ylabel('高通濾波器')
axis([0 0.5 -110 10])

% 將高斯白噪通過兩個濾波器
r=randn(16384,1);

x1=filter(b1,1,r);
x=x1(50:50+N);
[xpsd,F]=pwelch(x,N/4,1);
xpsd=10*log10(xpsd);
subplot(323);
plot(F,xpsd);
grid;
ylabel('x的功率譜')
axis([0 3 -70 10])

y1=filter(b2,1,r);
y=y1(50:50+N);
[ypsd,F]=pwelch(y,N/4,1);
ypsd=10*log10(ypsd);
subplot(324);
plot(F,ypsd);
grid;
ylabel('y的功率譜')
axis([0 3 -70 10])

% 估計 x和y 的相干函數；
[cxy,w]=cohere(x,y,N/4,Fs,hamming(N/4),0,'mean');
subplot(325);
plot(w,cxy);
grid;
ylabel('x與y的相干函數')

% 估計 x和y 的互功率譜；
[pxy,w]=csd(x,y,N/4,Fs,hamming(N/4),0,'mean');
pxy=20*log10(abs(pxy));
subplot(326);
plot(w,pxy);
grid;
ylabel('x與y的互功率譜')
axis([0 0.5 -70 10])