function [Xsig, fsig, X, f, E, Esig] = dft_signoi(x, T, Tm, fc, fmin, fmax)
% 071127 revised
[n, nch] = size(x);
X = fft(x) / n * 2;			% DFT
f = [0 : n-1]' / T;			% 周波数
m = round(T / Tm);			% 変調回数
ista = mod(find(f==fc),m);
isig = [ista : m : n]';		% 信号チャンネルの番号
Xsig = X(isig, :);			% 信号チャンネルのデータ
fsig = f(isig);				% 信号チャンネルの周波数
if nargin==6
    idx = fsig>=fmin & fsig<=fmax;	% 周波数の範囲を限定する場合
    Xsig = Xsig(idx, :);
    fsig = fsig(idx);
end

%%%% 一様なノイズレベル E を計算 %%%%
inoi = [1:n]';
if nargin==6
    inoi = inoi(f(inoi)>=fmin & f(inoi)<=fmax);	% 周波数の範囲を限定する場合
end
inoi = inoi(mod(inoi,m)~=ista);
E = sqrt(mean(abs(X(inoi,:)).^2, 1)/2);

%%%% 信号周波数ごとのノイズレベル Esig を計算 %%%%
isig = isig(f(isig)>=fmin & f(isig)<=fmax);
M = 5;
mm = floor(M/(m-1));
inoi = [-M-mm : M+mm];
inoi = inoi(mod(inoi,m)~=0);
inoi1 = repmat(isig,[1 size(inoi,2) nch]);
inoi2 = repmat(inoi,[size(isig,1) 1 nch]);
inoi3 = repmat(reshape([0:n:n*(nch-1)],[1 1 nch]),[size(isig,1) size(inoi,2) 1]);
inoi = inoi1 + inoi2;
inoi(inoi>n) = inoi(inoi>n) - n;
inoi(inoi<=0) = inoi(inoi<=0) + n;
inoi = inoi + inoi3;
Esig = squeeze(sqrt(mean(abs(X(inoi)).^2, 2)/2));
