%------begin computing the Least Square Estimate of intensity coefficients
frequency=[20;200;2000;20000];
timeIntvl=0:0.025:24.975;
timeIntvl=timeIntvl*(10^-3);  % transform the ms unit of time into s unit 
ctrlVarMatrix=sin(2*pi*(frequency*timeIntvl));   % compute the matrix of the control variables
fid=fopen('observationsOfSample.txt','r');
obstnVector_Y=fscanf(fid,'%f'); % read observations of sample from file
fclose(fid);
obstnVector_Y=obstnVector_Y*(10^-3);  % transform the mV unit of voice into V unit
factorLSE1=ctrlVarMatrix*transpose(ctrlVarMatrix);
if(rank(factorLSE1)<min(size(factorLSE1)))
    factorLSE2=inv(factorLSE1+eye(4,4)*(10^-7));  % modify the factorLSE1:add up factorLSE1 and a unit matrix that multiled by a tiny value  
else
    factorLSE2=inv(factorLSE1);  % compute the generalized inverse matrix of factorLSE1 
end
factorLSE=factorLSE2*ctrlVarMatrix*obstnVector_Y;  % compute the Least Square Estimate of intensity coefficients   
%------end computing the Least Square Estimate of intensity coefficients
%------begin computing the Least Square Estimate of variance
Qe=0;
for i=1:size(obstnVector_Y,1),
    Qe=Qe+(obstnVector_Y(i)-(transpose(ctrlVarMatrix(:,i))*factorLSE))^2;
end
varianceLSE=Qe/(size(obstnVector_Y,1)-size(frequency,1)-1);
%------end computing the Least Square Estimate of variance
%------output the result of Least Square Estimate
disp(['the LSE of factors from 1 to ',num2str(size(factorLSE,1)),' as follow : '])
disp(vpa(factorLSE,4))
disp('the LSE of variance as follow : ')
disp(vpa(varianceLSE,4))
plot(timeIntvl,obstnVector_Y)