0001
0002
0003
0004
0005
0006
0007
0008
0009 function varargout = mlpsd_m(varargin)
0010
0011 import utils.const.*
0012
0013 utils.helper.msg(msg.PROC1, 'using MATLAB to compute core DFT');
0014
0015
0016 x = varargin{1};
0017 f = varargin{2};
0018 r = varargin{3};
0019 m = varargin{4};
0020 L = varargin{5};
0021 fs = varargin{6};
0022 win = varargin{7};
0023 order = varargin{8};
0024 olap = varargin{9};
0025
0026 twopi = 2.0*pi;
0027
0028 nx = length(x);
0029 nf = length(f);
0030 Sxx = zeros(nf,1);
0031 S = zeros(nf,1);
0032 ENBW = zeros(nf,1);
0033
0034 disp_each = round(nf/100)*10;
0035
0036 for j=1:nf
0037
0038 if mod(j,disp_each) == 0 || j == 1 || j == nf
0039 utils.helper.msg(msg.PROC1, 'computing frequency %04d of %d: %f Hz', j, nf, f(j));
0040 end
0041
0042
0043 l = L(j);
0044 switch win.type
0045 case 'Kaiser'
0046 win = specwin(win.type, l, win.psll);
0047 otherwise
0048 win = specwin(win.type, l);
0049 end
0050
0051 p = 1i * twopi * m(j)/l.*[0:l-1];
0052 C = win.win .* exp(p);
0053 Xolap = (1-olap);
0054
0055 A = 0.0;
0056
0057
0058 segLen = l;
0059 nData = length(x);
0060 ovfact = 1 / (1 - olap);
0061 davg = (((nData - segLen)) * ovfact) / segLen + 1;
0062 navg = round(davg);
0063
0064
0065 if navg == 1
0066 shift = 1;
0067 else
0068 shift = (nData - segLen) / (navg - 1);
0069 end
0070 if shift < 1 || isnan(shift)
0071 shift = 1;
0072 end
0073
0074
0075
0076 start = 1.0;
0077 for ii = 1:navg
0078
0079 istart = round (start);
0080 start = start + shift;
0081
0082
0083 xs = [x(istart:istart+l-1)].';
0084
0085
0086 switch order
0087 case -1
0088
0089 case 0
0090 xs = xs - mean(xs);
0091 case 1
0092 xs = detrend(xs);
0093 otherwise
0094 xs = polydetrend(xs, order);
0095 end
0096
0097
0098 a = C*xs;
0099 A = A + a*conj(a);
0100
0101 end
0102
0103 if mod(j,disp_each) == 0 || j == 1 || j == nf
0104 utils.helper.msg(msg.PROC2, 'averaged %d segments', navg);
0105 end
0106
0107 A2ns = 2.0*A/navg;
0108 S1 = win.ws;
0109 S12 = S1*S1;
0110 S2 = win.ws2;
0111 ENBW(j) = fs*S2/S12;
0112 Sxx(j) = A2ns/fs/S2;
0113 S(j) = A2ns/S12;
0114
0115 end
0116
0117 varargout{1} = S;
0118 varargout{2} = Sxx;
0119 varargout{3} = ENBW;
0120
0121 end
0122