TSplot/PspecCB.m at main · magnetotellurics/TSplot · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
if(UpdateSpecPlot)
   action =  'SpecPlot';
   tmpfig = gcf;
   ii = TSfigInd(tmpfig.Number);
   if length(TSp) >= ii
      if ~isempty(TSp{ii})
         h = findobj;
         if(isempty(find(TSp{ii}.hAx==h)))
            UpdateSpecPlot = 0;
            return
         end
      else
         return
      end
   else
      return
   end
else
   action = get(gcbo,'Tag');
   ii = get(gcf,'UserData');
   h = gcbo;
end

switch action
   case 'SpecCheckBox'

   case 'SpecPlot'

      FullTS = get(findobj('Tag','FullTS','Parent',TSp{ii}.hFig),'value');
      Nwin = get(findobj('Tag','Nwindows','Parent',TSp{ii}.hFig),'string');
      Nwin = str2num(Nwin);
      hCh = findobj('Tag','SpecCheckBox','Parent',TSp{ii}.hFig);
      specCh = [];
      for k = length(hCh):-1:1
         if get(hCh(k),'value')
            l = get(hCh(k),'userdata');
            specCh = [ specCh l];
         end
      end
      ctitle1 = ['Power Spectra'];
      lgnd = [];
      lineS = {'b-','r-','g-','k-','m-','b--','r--','g--','k--','m--'};
      hLine = zeros(length(specCh),1);
      for k = 1:length(specCh)
         ich = specCh(k);
         lgnd{k} = [ TSd{ii}.ch_id(ich,:) TSd{ii}.sta(ich,:)];
         if FullTS
            d = TSdata{ii}(ich,:) ; %./TSd1.scale(ich);  % gain corr.
         else
            i0 = TSw{ii}.i0;
            i1 = TSw{ii}.i1;
            d = TSdata{ii}(ich,i0:i1) ; %./TSd1.scale(ich); % gain corr.
         end
         %   prune off NaN's at start begining
         isnd = isnan(d);
         if any(isnd)
            iGood = find(1-isnd);
            dBar = mean(d(iGood));
            d = d - dBar;
            clear iGood;
            d(isnd) = 0;
         else
            d = d - mean(d);
         end
         clear isnd;
         d = quickClean(d,2,10,.1);

         nData = length(d);
         window = floor(nData/Nwin);
         %  make sure window isn't too weird for fft
         %    (this will cause part of time series to be discarded)
         window = 8*floor(window/8);
%         [ps,w] = pwelch(d,window);
         nfft = window  ;             % TH psd
         noverlap = nfft/2;
         fs = 1/(TSd{ii}.dt);
%        w = frequency in Hz
         [ps,w]= pwelch(d,hanning(nfft),noverlap,nfft,fs,'onesided');
%         w = w/(TSd{ii}.dt*2*pi);
         ps = sqrt(ps);
         hax = TSp{ii}.hAx;
         axes(hAx);
         hLine(k) = loglog(w(2:length(w)),ps(2:length(ps)),lineS{k});
         hold on
         set(hLine(k),'LineWidth',2);
      end
      if length(specCh) > 0
         set(hax,'FontWeight','demi');
         xlabel('Hz')
%         ylabel('Amplitude (nT/Hz^{1/2} or mV/Hz^{1/2})');
         ylabel('Amplitude (V/Hz^{1/2})');
         title(ctitle1);
         legH = legend(hLine,lgnd,'Location','Best');
         set (legH,'Interpreter','none')
         hold off
      end

   case 'NWin'
      NWin = str2num(get(h,'String'));
      if NWin > 0
         TSp{ii}.NWin=NWin;
      else
         set(h,'String',num2str(TSp{ii}.NWin));
      end
end
UpdateSpecPlot = 0;