急！！下面的是一个时频分析的wigner-ville分布的matlab代码，求高手转换成C++，我用的是vs2003做MFC项目

[cjohn1988]

如下代码，其中的tftb_window窗采用h=0.54 - 0.46*cos(2.0*pi*(1:N)'/(N+1));
其他的函数应该基本上没有大的影响的
function [tfr,t,f] = tfrsp(x,t,N,h,trace);
%TFRSP Spectrogram time-frequency distribution.
% [TFR,T,F]=TFRSP(X,T,N,H,TRACE) computes the Spectrogram
% distribution of a discrete-time signal X.
%
% X : signal.
% T : time instant(s) (default : 1:length(X)).
% N : number of frequency bins (default : length(X)).
% H : analysis window, H being normalized so as to
% be of unit energy. (default : Hamming(N/4)).
% TRACE : if nonzero, the progression of the algorithm is shown
% (default : 0).
% TFR : time-frequency representation. When called without
% output arguments, TFRSP runs TFRQVIEW.
% F : vector of normalized frequencies.
%
% Example :
% sig=fmlin(128,0.1,0.4);
% h=tftb_window(17,'Kaiser'); tfrsp(sig,1:128,64,h,1);
%
% [tfr,t,freq]=tfrsp(sig,1:128,64,h,1); plot(fftshift(freq),fftshift(tfr(:,100)))
%
% See also all the time-frequency representations listed in
% the file CONTENTS (TFR*)

if (nargin == 0),
error('At least 1 parameter required');
end;
[xrow,xcol] = size(x);
if (nargin < 1),
error('At least 1 parameter is required');
elseif (nargin <= 2),
N=xrow;
end;

hlength=floor(N/4);
hlength=hlength+1-rem(hlength,2);

if (nargin == 1),
t=1:xrow; h = tftb_window(hlength); trace=0;
elseif (nargin == 2) | (nargin == 3),
h = tftb_window(hlength); trace=0;
elseif (nargin == 4),
trace = 0;
end;

if (N<0),
error('N must be greater than zero');
end;
[trow,tcol] = size(t);
if (xcol==0) | (xcol>2),
error('X must have one or two columns');
elseif (trow~=1),
error('T must only have one row');
elseif (2^nextpow2(N)~=N),
fprintf('For a faster computation, N should be a power of two\n');
end;

[hrow,hcol]=size(h); Lh=(hrow-1)/2;
if (hcol~=1)|(rem(hrow,2)==0),
error('H must be a smoothing window with odd length');
end;

tfr= zeros (N,tcol) ;
if trace, disp('Spectrogram'); end;
for icol=1:tcol,
ti= t(icol); tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
indices= rem(N+tau,N)+1;
if trace, disprog(icol,tcol,10); end;
tfr(indices,icol)=x(ti+tau).*conj(h(Lh+1+tau))/norm(h(Lh+1+tau));
end;
tfr=abs(fft(tfr)).^2;
if trace, fprintf('\n'); end;

if (nargout==0),
tfrqview(tfr,x,t,'tfrsp',h);
elseif (nargout==3),
if rem(N,2)==0,
f=[0:N/2-1 -N/2:-1]'/N;
else
f=[0:(N-1)/2 -(N-1)/2:-1]'/N;
end;
end;

[cjohn1988]

求各位大牛帮帮忙吧~~