%文件名:csf.m
%程序員:郭遲
%編寫時間:2004.3.22
%函數功能:本函數將完成一個CSF的設計
%輸入格式舉例:filtercoefficients=csf
%參數說明:
%filtercoefficients為CSF的濾波器系數
function filtercoefficients=csf()

%調用子函數計算頻率響應矩陣
Fmatrix = csfmat;
%畫出頻率響應
%figure(1);mesh(Fmatrix),title('頻率響應'),xlabel('水平方向空間頻率');ylabel('垂直方向空間頻率');zlabel('CSF頻率響應'); 
%利用FSAMP2函數計算頻率系數
filtercoefficients=fsamp2(Fmatrix);   

%===================
%子函數,計算頻率響應矩陣
function Fmatrix=csfmat()
%===================
u = -20:1:20; 
v = -20:1:20;
n = length(u);
Z = zeros(n);
for i=1:n
	for j=1:n
		Z(i,j)=csffun(u(i),v(j));	% 調用子函數計算響應空間頻率下的頻響
	end
end
Fmatrix = Z;

%========================
%子函數,計算u,v下的頻率響應
function Sa=csffun(u,v)
%========================
%CSF頻率響應
sigma = 2;
f = sqrt(u.*u+v.*v);
w = 2*pi*f/60;
Sw = 1.5*exp(-sigma^2*w^2/2)-exp(-2*sigma^2*w^2/2);

%高頻修正
sita = atan(v./(u+eps));% eps= 2^(-52)是避免0的一種修正
bita = 8;
f0 = 11.13;
w0 = 2*pi*f0/60;
Ow = ( 1 + exp(bita*(w-w0)) * (cos(2*sita))^4) / (1+exp(bita*(w-w0)));

%最終結果
Sa = Sw * Ow;