Contents

function g2(t,y,P)

% Starting from a certain signal that you created using g1, attempt to
% recover the amplitudes and periods from the composite signal
%
% How to run:
%

OPTION 1:

Generate a signal of your choice, e.g.

[t,y]=g1([1 4 5],[10 30 40],[0 pi/2 0]);

Attempt to recover the amplitudes GIVEN the periods - put in

"reasonable guesses'

g2(t,y,[ 9 20 42]); g2(t,y,[10 30 40]); g2(t,y,[10 30 46]); % <------ A

The output is to screen, and to a plot

OPTION 2: NOISY DATA

Add some noise before you try this!

yn=y+randn(size(y))*std(y)/sqrt(5); g2(t,yn,[ 9 20 42]); g2(t,yn,[10 30 40]); % <------ B g2(t,yn,[10 30 46]);

COMPARE A with B - slightly off in period vs slighly noisy

OPTION 3: MISSING DATA, SET TO ZERO

Knock out some values at random

frax=3; ko=max(1,round(rand(round(length(y)/frax),1)*length(y))); yk=y; yk(ko)=0; g2(t,yk,[ 9 20 42]); g2(t,yk,[10 30 40]); % <------ C g2(t,yk,[10 30 46]);

OPTION 4: MISSING DATA, SIMPLY OMITTED

allt=1:length(t); allt(ko)=NaN; allt=allt(~isnan(allt)); g2(t(allt),y(allt),[ 9 20 42]); g2(t(allt),y(allt),[10 30 40]); % <------ D g2(t(allt),y(allt),[10 30 46]);

OPTION 5: INTERPOLATE MISSING DATA

yi=interp1(t(allt),y(allt),t,'spline'); 'linear' etc

Take out any NaNs that are interpolation artifacts

ti=t(~isnan(yi)); yi=yi(~isnan(yi)); g2(ti,yi,[ 9 20 42]); g2(ti,yi,[10 30 40]); % <------ E g2(ti,yi,[10 30 46]);

% Set your netid lab name
netid='mynetidl05';

% First of all, make sure that all inputs are column vectors
t=t(:);
y=y(:);

% The "model" is that the signal IS given by a weighted set of sines and
% cosines for whom you have just guessed the periods, but you don't know
% about the amplitudes. So this here is the 'design' matrix
% We are doing this all at the same time, folks, so check out what the
% next line does. 1./P are the frequencies. Do it for ONE period and you
% should understand it right away. Do it for TWO OR MORE and you'll love it
argm=[1./P(:)*t(:)']';
% This here is how you think the y-vector relates to the t-vector
% The first column is to allow for an offset - a nonzero mean
F=[ones(size(y)) cos(2*pi*argm) sin(2*pi*argm)];
% If you are unwilling to phrase the phasing this way, use complex
% exponentials and deal with the complex numbers that come out.

% And this here is how you "invert" the relationship to come up with the
% weighting constants that are the amplitudes for each of the terms, and by
% 'applying' the inverse to the data vector - a 'linear regression'!
A=pinv(F)*y;

% If this all were to be exact, the PREDICTION would be spot on (we'll see)
ypred=F*A;

% Calculate the residuals
resd=y-ypred;

% Calculate the root-mean-squared error
rmser=sqrt(sum(resd.^2)/length(resd));

% Compare the variance of the residual to the variance of the data: if
% the fit does a good job, the variance of the residuals should be a
% small fraction of the variance of the data
vard=100*[1-var(resd)/var(y)];

% Now remember that you need to rearrange the amplitudes to turn them
% into magnitudes that belong to a certain period for either the sine or
% the cosine, regardless of the phase angle. Look at the above to realize
% that we need to sum them in pairs starting from the second one

% Pull out the offset
m=A(1);
% Rearrange the amplitudes pairwise, sum, square, square-root
A=sqrt(sum(reshape(A(2:end),length(P),[]).^2,2));

% Report to screen, repeat the periods which you yourself put in
ostr='%5.2f  ';
fstr=repmat(ostr,1,length(A));
disp(sprintf(['\n amplitudes A ' fstr '\n periods    P ' fstr],...
             A,P))
disp(sprintf([' offset term  ' ostr],m))
disp(sprintf([' variance rd  ' deblank(ostr) '%%'],vard))

% Plot what you think is going on
figure
% Top panel %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(211)
% Plot the signal
plot(t,y,'k','linew',2)
% Plot the prediction of the data
hold on
plot(t,ypred,'b','linew',1)
hold off
% Embellish
grid on
xlabel('time [s]')
ylabel('signal and its prediction')
title(sprintf(['rms prediction error is %5.2f data units\n' ...
               'variance reduction ' deblank(ostr) '%%'],...
              rmser,vard))

% Bottom panel %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(212)
stem(P,A,'filled')
% Again, open up the axes a bit
ylim([-0. 1]*max(A)*1.1)
xlim([0 1]*max(P)*1.1)
% Embellish
grid on
xlabel('period [s]')
ylabel('amplitude [units of the data]')
title(sprintf(['%s - A ' fstr 'P ' fstr],...
               netid,A,P))

Not enough input arguments.

Error in g2 (line 61)
t=t(:);

Published with MATLAB® R2015b