% main parameters
 Fs = 4; % sampling frequency in case of
 uniform sampling
 startTime = 1; % start of sampling time
 endTime = 400; % end of sampling time
 % Note 1/100 Hz resolution will be
% obtained only if the time period
 is
% longer than 100 seconeds.
 % frequencies of y(t) and x(t)
 yFrequencies = [1.01,1.02];
 yPhases = [0 ,0 ];
 xFrequencies = 0.05:0.01:2.00; % you may try other frequencies
 %xPhases = 2*pi*rand(size(xFrequencies)); % you may try this
 xPhases = zeros(size(xFrequencies));
%--------------------------------------------------------------------------
Building samples
-------------------------------------------------------
 if Ctrl.DoUniformSampling == true
2
 % build time vector and find number of samples
 t = (startTime:1/Fs:endTime)'; % time vector
 N = length(t); % number of samples
 else
 % creat non-uniform samples
 N = Fs*(endTime-startTime); % number of samples
 t = startTime + rand(N,1)*(endTime-startTime);
 t = sort(t);
 end
 % build samples of y(t)
 yn = zeros(N,1); % initialize with zer
 for i = 1 : length(yFrequencies)
 yn = yn + cos(2*pi*yFrequencies(i).*t+yPhases(i));
 end
 % build samples of x(t)
 xn = zeros(N,1); % initialize with zer
 for i = 1 : length(xFrequencies)
 xn = xn + cos(2*pi*xFrequencies(i).*t+xPhases(i));
 end
 % build samples of y' (these are measurements)
 ypn = yn + xn;
%--------------------------------------------------------------------------
Linear least sqaure estimation
-----------------------------------------
For linear estimation we have to expany the cos function in y''. In this case we have cos(a+b) = cos(a)cos(b)
- sin(a)sin(b). Therefor, y'' = A*cos(2*pi*f*t+p) => y'' = A*cos(2*pi*f*t)*cos(p) -A*sin(2*pi*f*t)*sin(p)
=> y'' = A*cos(p)*cos(2*pi*f*t) -A*sin(p)*sin(2*pi*f*t) => y'' = B *cos(2*pi*f*t) -C *sin(2*pi*f*t)
where B = A*cos(p) and C = A*sin(p)
 % now we are going for all frequecies in x estimate the amplitude
 and
 % phase from the samples of y'
 % In matrix form for each row we have
 % ynp[n] = [cos(2*pi*f*t[n]) -sin(2*pi*f*t[n])] * [B C]'
 % and for all rows we have
 % Yp(N x 1) = H(N x 2) * transpose([B C])
 % where H is a matrix built from [cos(2*pi*f*t[n]) -
sin(2*pi*f*t[n])]
 % for all sampling times
 % define two vectors to store the results
 linearEstAmps = zeros(length(xFrequencies),1);
 linearEstPhases = zeros(length(xFrequencies),1);
 for i = 1 : length(xFrequencies)
3
 % Yp is already calculated
 Yp = ypn;
 % building H
 H = [cos(2*pi*xFrequencies(i)*t) -
sin(2*pi*xFrequencies(i)*t)];
 % least square estimation
 % ref: https://en.wikipedia.org/wiki/
Linear_least_squares_(mathematics)#Motivational_example
 %estBC = (transpose(H) * H)^-1*transpose(H) * Yp;
 estBC = pinv(H) * Yp; % due to numerical error it is
 better to use pinv
 % find B and C estimate
 estB = estBC(1);
 estC = estBC(2);
 % find amplitude and phase
 estA = sqrt(estB.^2 + estC.^2); % as A = sqrt(B^2 + C^2)
 estP = atan2(estC,estB); % as P = atan(C/B)
 % store the results
 linearEstAmps(i) = estA;
 linearEstPhases(i) = estP;
 end
 % plot amplitude and phase
 figure(1)
 clf
 subplot(2,1,1)
 bar(xFrequencies,linearEstAmps,'b')
 xlabel('frequencies')
 ylabel('estimated amplitude')
 grid on
 grid minor
 xlim([xFrequencies(1) xFrequencies(end)])
 hold on
 title('performance of linear LS (expect is 1)')
 subplot(2,1,2)
 bar(xFrequencies,linearEstPhases,'b')
 xlabel('frequencies')
 ylabel('estimated phase [rad]')
 title('(expect is 0)')
 grid on
 grid minor
 xlim([xFrequencies(1) xFrequencies(end)])
%--------------------------------------------------------------------------
4
Non-Linear least sqaure estimation
-------------------------------------
as we know y'' is a non-linear function of the phase parameter. Therefore, we need to use non-linear LS
algorithm, and solve the problem iteratively. we need to find the derivaties of y'' w.r.t. the A and P dy''/
dA = cos(2*pi*f*t+P) dy''/dP = -A*sin(2*pi*f*t+P)
 numOfIteration = 70;
 convergenceRho = 0.2;
 % define two vectors to store the results
 nonLinearEstAmps = zeros(length(xFrequencies),1);
 nonLinearEstPhases = zeros(length(xFrequencies),1);
 for i = 1 : length(xFrequencies)
 % Yp is already calculated
 Yp = ypn;
 % Initialize estimation point
 A = 1;
 P = 0;
 for k = 1 : numOfIteration
 % build y'' samples from the estimation
 Yppn = A*cos(2*pi*xFrequencies(i)*t+P);
5
 % find residual (measurements - estimated samples)
 r = Yp - Yppn;
 % build Jacobian matrix J = [df/dx1 df/dx2]
 % reference: https://en.wikipedia.org/wiki/
Jacobian_matrix_and_determinant
 J = [cos(2*pi*xFrequencies(i)*t+P),-
A*sin(2*pi*xFrequencies(i)*t+P)];
 % find new A and P
 estAP = [A;P] + convergenceRho.*pinv(J)*r;
 % update A and P
 A = estAP(1);
 P = estAP(2);
 % remove ambiguities
 if A < 0
 A = -A;
 P = P + pi;
 end
 end
 % store the results
 nonLinearEstAmps(i) = A;
 nonLinearEstPhases(i) = P;
 end
 % plot amplitude and phase
 figure(2)
 clf
 subplot(2,1,1)
 bar(xFrequencies,nonLinearEstAmps,'b')
 xlabel('frequencies')
 ylabel('estimated amplitude')
 grid on
 grid minor
 xlim([xFrequencies(1) xFrequencies(end)])
 hold on
 title('performance of non-linear LS using Gauss-Newton')
 subplot(2,1,2)
 bar(xFrequencies,nonLinearEstPhases,'b')
 xlabel('frequencies')
 ylabel('estimated phase [rad]')
 grid on
 grid minor
 xlim([xFrequencies(1) xFrequencies(end)])
6
Analysis
---------------------------------------------------------------
find estimate of y' from estimated frequencies and phases
 estLinearYp = zeros(N,1);
 estNonLinearYp = zeros(N,1);
 for i = 1 : length(xFrequencies)
 estLinearYp = estLinearYp + ...
 linearEstAmps(i)*cos(2*pi*xFrequencies(i)*t +
 linearEstPhases(i));
 estNonLinearYp = estNonLinearYp + ...
 nonLinearEstAmps(i)*cos(2*pi*xFrequencies(i)*t +
 nonLinearEstPhases(i));
 end
 figure(3)
 plot(t,ypn,'b',...
 t,estLinearYp,'r--',...
 t,estNonLinearYp,'k-.')
 xlabel('time in sec')
 ylabel('amplitude')
 grid on
 grid minor
 title('performance of the estimation')
7
 legend('Actual','linear-est','non-linear-est')
 % find sum error
 estLinearError_dB = 10*log10(sum(abs(ypn - estLinearYp
 ).^2));
 estNonLinearError_dB = 10*log10(sum(abs(ypn -
 estNonLinearYp).^2));
 fprintf('Estimation error: n tLinear = %0.2f dBntNonlinear
= %0.2f dBn',...
 estLinearError_dB,estNonLinearError_dB)
Estimation error:
 Linear = 12.21 dB
 Non-linear = 12.21 dB