home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / CS2.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  5.0 KB  |  159 lines
  1. echo on,clc
  2. % CASE STUDY # 2 ANALYSING A SIGNAL FROM A TRANSFORMER
  3. %
  4. % In this case study we shall consider the current signal from the
  5. % R-phase when a  400 kV three-phase transformer is energized.
  6. % The measurements were performed by Sydkraft AB in Sweden.
  7.  
  8. load current.mat
  9.  
  10. % The sampling interval is 1 ms.
  11. % Let's take a look at the  data.
  12.  
  13. pause,idplot(i4r,[],0.001),pause
  14.  
  15. % A close up:
  16.  
  17. pause,idplot(i4r,201:250,0.001),pause
  18.  
  19. % Let us first compute the periodogram:
  20.  
  21. ge = etfe(i4r);
  22. ge=sett(ge,0.001); % Adjusting to the sampling interval 1 ms
  23.  
  24. pause,bodeplot(ge),pause
  25.  
  26. % A smoothed periodogram is obtained by
  27.  
  28. ges=etfe(i4r,length(i4r)/4);ges=sett(ges,0.001);
  29.  
  30. % A plot with linear frequency scale  is given by
  31.  
  32. pause, ffplot(ges),grid, pause
  33.  
  34. % We clearly see the dominant frequency component of 50 Hz, and its harmonics.
  35.  
  36. % Comparing with the pure periodogram gives:
  37.  
  38. pause, ffplot([ges ge]),pause
  39.  
  40. % The standard spectral analysis estimate gives (with the default window 
  41. % size, which is not adjusted to resonant spectra)
  42.  
  43. gs =  spa(i4r);
  44. gs = sett(gs,0.001); % Adjusting to the sampling interval 1 ms.
  45. pause,ffplot([gs ges]),pause
  46.  
  47. % We see that a very large lag window will be required to see all the fine
  48. % resonances of the signal. Standard spectral analysis does not work well.
  49.  
  50. % Let us instead compute spectra by parametric AR-methods. Models of 2nd
  51. % 4th and 8th order are obtained by
  52.  
  53. t2=ar(i4r,2);t2=sett(t2,0.001); % (The sampling interval)
  54. t4=ar(i4r,4);t4=sett(t4,0.001);
  55. t8=ar(i4r,8);t8=sett(t8,0.001);
  56. g2=th2ff(t2);g4=th2ff(t4);g8=th2ff(t8);
  57.  
  58. % Let us take a look at the spectra:
  59.  
  60. pause,ffplot([g2 g4 g8 ges]),pause
  61.  
  62. % We see that the parametric spectra are not capable of picking up the
  63. % harmonics. The reason is that the AR-models attach too much attention to
  64. % the higher frequencies, which are difficult to  model. (See Ljung (1987)
  65. % Example 8.5)
  66. % We will have to go to high order models before the harmonics are picked
  67. % up.
  68.  
  69. % What will a useful order be?
  70.  
  71. V=arxstruc(i4r(1:301),i4r(302:601),[1:30]'); % Checking all order up to 30
  72. pause,nn=selstruc(V,'log');
  73.  
  74. % We see a dramatic drop for n=20, so let's pick that order
  75.  
  76. t20=ar(i4r,20);t20=sett(t20,0.001);g20=trf(t20);
  77.  
  78. pause,ffplot([ges g20]),pause
  79.  
  80. % All the harmonics are now picked up, but why has the level dropped?
  81. % The reason is that t20 contains very thin but high peaks. With the
  82. % crude gitter of freqeuncy points in g20 we simply don't see the 
  83. % true levels of the peaks. We can illustrate this as follows:
  84.  
  85. g20c=trf(t20,[],[551:550+128]/600*150*2*pi); % A frequency region around
  86.                                              % 150 Hz
  87.  
  88. pause,ffplot([ges g20 g20c]),pause
  89.  
  90. % If we are primarily interested in the lower harmonics, and want to
  91. % use lower order models we will have to apply prefiltering of
  92. % the data. We select a fifth order Butterworh filter with cut-off
  93. % frequency at 200 Hz. (This should cover the 50, 100 and 150 Hz modes):
  94.  
  95. [bb,ba]=butter(5,200/500); % 500 Hz is the Nyquist frequency
  96. i4rf = filtfilt(bb,ba,i4r);
  97.  
  98. t8f=ar(i4rf,8); t8f=sett(t8f,0.001);g8f=th2ff(t8f);
  99.  
  100. % Let us now compare the spectrum obtained from the filtered data (8th
  101. % order model) with that for unfiltered data (8th order) and with the
  102. % periodogram:
  103.  
  104. pause,ffplot([g8f g8 ges]),pause
  105.  
  106. % We see that with the filtered data we pick up the three first peaks in
  107. % the spectrum quite well. 
  108.  
  109. % We can compute the numerical values of the resonances as follows:
  110. % The roots of a sampled sinusoid of frequency om are located on
  111. % the unit circle at exp(i*om*T), T being the sampling interval. We
  112. % thus proceed as follows:
  113.  
  114. pause,a=th2poly(t8f) % The AR-polynomial
  115. omT=angle(roots(a))'
  116. freqs=omT/0.001/2/pi'  % In Hz
  117. pause
  118.  
  119. % We thus find the harmonics quite well.  We could also test how well
  120. % the model t8f is capable of predicting the signal, say 50 ms ahead:
  121.  
  122. pause,compare(i4rf,t8f,50,201:500);pause
  123.  
  124. % If we were interested in the fourth and fifth harmonics (around
  125. % 200 and 250 Hz) we would proceed as follows:
  126.  
  127. [bb,ba]=butter(5,[150 300]/500);
  128. i4rff=filtfilt(bb,ba,i4r);
  129. t8f=ar(i4rff,8); t8f=sett(t8f,0.001);
  130. g8f=th2ff(t8f);
  131.  
  132. pause,ffplot([ges g8f]),pause
  133.  
  134. % We thus see that with proper prefiltering, low order parametric models
  135. % can be built that give good descriptions of the signal over desired
  136. % frequency ranges.
  137.  
  138. % What would the 20th order model give in terms of estimated harmonics:
  139.  
  140. pause,a=th2poly(t20) % The AR-polynomial
  141. omT=angle(roots(a))'
  142. freqs=omT/0.001/2/pi'  % In Hz
  143. pause
  144.  
  145. % We see that this model picks  up the harmonics very well.
  146. % This model will also predict as follows:
  147.  
  148. pause,compare(i4r,t20,50,201:500);pause
  149.  
  150. % For a complete model of the signal, t20 thus is the natural choice,
  151. % both in terms of finding the harmonics and in prediction capabilities.
  152. % For models in certain frequency ranges we can however do very well
  153. % with lower order models, but we then have to prefilter the data
  154. % accordingly.
  155.  
  156. pause
  157. echo off
  158.  
  159.