0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 function varargout = mlpsd_mex(varargin)
0011
0012 import utils.const.*
0013
0014 utils.helper.msg(msg.PROC1, 'using ltpda_dft.mex to compute core DFT');
0015
0016
0017 x = varargin{1};
0018 f = varargin{2};
0019 r = varargin{3};
0020 m = varargin{4};
0021 L = varargin{5};
0022 fs = varargin{6};
0023 win = varargin{7};
0024 order = varargin{8};
0025 olap = varargin{9};
0026 Lmin = varargin{10};
0027
0028 twopi = 2.0*pi;
0029
0030 nx = length(x);
0031 nf = length(f);
0032 Sxx = zeros(nf,1);
0033 S = zeros(nf,1);
0034 ENBW = zeros(nf,1);
0035
0036 disp_each = round(nf/100)*10;
0037
0038 minReached = 0;
0039
0040 for j=1:nf
0041
0042 if mod(j,disp_each) == 0 || j == 1 || j == nf
0043 utils.helper.msg(msg.PROC1, 'computing frequency %04d of %d: %f Hz', j, nf, f(j));
0044 end
0045
0046
0047 l = L(j);
0048
0049
0050
0051 if ~minReached
0052 switch win.type
0053 case 'Kaiser'
0054 win = specwin(win.type, l, win.psll);
0055 otherwise
0056 win = specwin(win.type, l);
0057 end
0058 if l==Lmin
0059 minReached = 1;
0060 end
0061 end
0062
0063 p = 1i * twopi * m(j)/l.*(0:l-1);
0064 C = win.win .* exp(p);
0065
0066
0067 [A, nsegs] = ltpda_dft(x, l, C, olap, order);
0068
0069 if mod(j,disp_each) == 0 || j == 1 || j == nf
0070 utils.helper.msg(msg.PROC2, 'averaged %d segments', nsegs);
0071 end
0072
0073 A2ns = 2.0*A/nsegs;
0074 S1 = win.ws;
0075 S12 = S1*S1;
0076 S2 = win.ws2;
0077 ENBW(j) = fs*S2/S12;
0078 Sxx(j) = A2ns/fs/S2;
0079 S(j) = A2ns/S12;
0080
0081 end
0082
0083 varargout{1} = S;
0084 varargout{2} = Sxx;
0085 varargout{3} = ENBW;
0086
0087 end
0088