function [fCFA, bCFA, T] = fs( sCFA )
% [fCFA, bCFA, T] = fs( sCFA )
%      to find the frequency structure for a given CFA, sCFA.
%
% sCFA : a matrix of a CFA pattern with R, G and B as symbolic variables.
% fCFA : a matrix of the same size and for the frequency structure of sCFA.
% bCFA : a boolean matrix to show the nonzero components in fCFA.
% T    : the multiplexing matrix for the frequency structure of sCFA.
%
%   For any commercial use of this program, a written permission can be 
%   obtained from the authors. please email: phao@eecs.qmul.ac.uk
%
%   If this program is used for research, please cite in your publications:
%
%   * Pengwei Hao, Yan Li, Zhouchen Lin, and Eric Dubois, "A Geometric
%     Method for Optimal Design of Color Filter Arrays", IEEE Transactions
%     on Image Processing, Vol. 20, No. 3, pp. 709-722, March 2011.
%
%   * Yan Li, Pengwei Hao, and Zhouchen Lin, "The Frequency Structure Matrix: 
%     A Representation of Color Filter Arrays", International Journal 
%     of Imaging Systems and Technology, Vol.21, No.1, pp.101-106, 2011.
%
%   By Pengwei Hao and Yan Li, 2007-2012

[m,n] = size(sCFA);

syms R G B real

% 2D symbolic Fourier transform
fCFA = simple( (fft(eye(m))/m) * sCFA * (fft(eye(n))/n).' );

% To set the small components to 0
th = 2^-16; % threshold to find the nonzero components
fR = subs(fCFA,{R,G,B},{1,0,0});
fG = subs(fCFA,{R,G,B},{0,1,0});
fB = subs(fCFA,{R,G,B},{0,0,1});
bCFA = double( abs(fR) + abs(fG) + abs(fB) ) > th;

% To show fCFA simpler and remove the small components
fCFA = simple( (R*fR + G*fG + B*fB) .* bCFA );

% To find the multiplexing matrix
[row,col] = find(bCFA);
T = [];
for k = 1:length(row),
    fnz = fCFA(row(k),col(k));
    T = [T; subs(fnz,{R,G,B},{1,0,0}), subs(fnz,{R,G,B},{0,1,0}), subs(fnz,{R,G,B},{0,0,1}) ];
end;

end
