echo on,clc
% CASE STUDY # 2 ANALYSING A SIGNAL FROM A TRANSFORMER
%
% In this case study we shall consider the current signal from the
% R-phase when a  400 kV three-phase transformer is energized.
% The measurements were performed by Sydkraft AB in Sweden.

load current.mat

% The sampling interval is 1 ms.
% Let's take a look at the  data.

pause,idplot(i4r,[],0.001),pause

% A close up:

pause,idplot(i4r,201:250,0.001),pause

% Let us first compute the periodogram:

ge = etfe(i4r);
ge=sett(ge,0.001); % Adjusting to the sampling interval 1 ms

pause,bodeplot(ge),pause

% A smoothed periodogram is obtained by

ges=etfe(i4r,length(i4r)/4);ges=sett(ges,0.001);

% A plot with linear frequency scale  is given by

pause, ffplot(ges),grid, pause

% We clearly see the dominant frequency component of 50 Hz, and its harmonics.

% Comparing with the pure periodogram gives:

pause, ffplot([ges ge]),pause

% The standard spectral analysis estimate gives (with the default window 
% size, which is not adjusted to resonant spectra)

gs =  spa(i4r);
gs = sett(gs,0.001); % Adjusting to the sampling interval 1 ms.
pause,ffplot([gs ges]),pause

% We see that a very large lag window will be required to see all the fine
% resonances of the signal. Standard spectral analysis does not work well.

% Let us instead compute spectra by parametric AR-methods. Models of 2nd
% 4th and 8th order are obtained by

t2=ar(i4r,2);t2=sett(t2,0.001); % (The sampling interval)
t4=ar(i4r,4);t4=sett(t4,0.001);
t8=ar(i4r,8);t8=sett(t8,0.001);
g2=th2ff(t2);g4=th2ff(t4);g8=th2ff(t8);

% Let us take a look at the spectra:

pause,ffplot([g2 g4 g8 ges]),pause

% We see that the parametric spectra are not capable of picking up the
% harmonics. The reason is that the AR-models attach too much attention to
% the higher frequencies, which are difficult to  model. (See Ljung (1987)
% Example 8.5)
% We will have to go to high order models before the harmonics are picked
% up.

% What will a useful order be?

V=arxstruc(i4r(1:301),i4r(302:601),[1:30]'); % Checking all order up to 30
pause,nn=selstruc(V,'log');

% We see a dramatic drop for n=20, so let's pick that order

t20=ar(i4r,20);t20=sett(t20,0.001);g20=trf(t20);

pause,ffplot([ges g20]),pause

% All the harmonics are now picked up, but why has the level dropped?
% The reason is that t20 contains very thin but high peaks. With the
% crude gitter of freqeuncy points in g20 we simply don't see the 
% true levels of the peaks. We can illustrate this as follows:

g20c=trf(t20,[],[551:550+128]/600*150*2*pi); % A frequency region around
                                             % 150 Hz

pause,ffplot([ges g20 g20c]),pause

% If we are primarily interested in the lower harmonics, and want to
% use lower order models we will have to apply prefiltering of
% the data. We select a fifth order Butterworh filter with cut-off
% frequency at 200 Hz. (This should cover the 50, 100 and 150 Hz modes):

[bb,ba]=butter(5,200/500); % 500 Hz is the Nyquist frequency
i4rf = filtfilt(bb,ba,i4r);

t8f=ar(i4rf,8); t8f=sett(t8f,0.001);g8f=th2ff(t8f);

% Let us now compare the spectrum obtained from the filtered data (8th
% order model) with that for unfiltered data (8th order) and with the
% periodogram:

pause,ffplot([g8f g8 ges]),pause

% We see that with the filtered data we pick up the three first peaks in
% the spectrum quite well. 

% We can compute the numerical values of the resonances as follows:
% The roots of a sampled sinusoid of frequency om are located on
% the unit circle at exp(i*om*T), T being the sampling interval. We
% thus proceed as follows:

pause,a=th2poly(t8f) % The AR-polynomial
omT=angle(roots(a))'
freqs=omT/0.001/2/pi'  % In Hz
pause

% We thus find the harmonics quite well.  We could also test how well
% the model t8f is capable of predicting the signal, say 50 ms ahead:

pause,compare(i4rf,t8f,50,201:500);pause

% If we were interested in the fourth and fifth harmonics (around
% 200 and 250 Hz) we would proceed as follows:

[bb,ba]=butter(5,[150 300]/500);
i4rff=filtfilt(bb,ba,i4r);
t8f=ar(i4rff,8); t8f=sett(t8f,0.001);
g8f=th2ff(t8f);

pause,ffplot([ges g8f]),pause

% We thus see that with proper prefiltering, low order parametric models
% can be built that give good descriptions of the signal over desired
% frequency ranges.

% What would the 20th order model give in terms of estimated harmonics:

pause,a=th2poly(t20) % The AR-polynomial
omT=angle(roots(a))'
freqs=omT/0.001/2/pi'  % In Hz
pause

% We see that this model picks  up the harmonics very well.
% This model will also predict as follows:

pause,compare(i4r,t20,50,201:500);pause

% For a complete model of the signal, t20 thus is the natural choice,
% both in terms of finding the harmonics and in prediction capabilities.
% For models in certain frequency ranges we can however do very well
% with lower order models, but we then have to prefilter the data
% accordingly.

pause
echo off

