diff --git a/mfiles/contwtgn.m b/mfiles/contwtgn.m deleted file mode 100644 index 566922e..0000000 --- a/mfiles/contwtgn.m +++ /dev/null @@ -1,117 +0,0 @@ -function [scalo,f,T,a,wt,wavescaled] = contwtgn(x,fmin,fmax,N,wave); - -% CONTWTGN [scalo,f,T,a,wt,wavescaled] = contwtgn(x,fmin,fmax,N,wave) -% computes a continuous wavelet transform. -% -% Input: -% x signal (in time) to be analyzed -% [fmin,fmax] respectively lower and upper frequency bounds of -% the analysis (in cycles/sec). -% [N] number of analyzed voices -% [wave] specifies the analyzing wavelet -% An order "wave" derivative of the Gaussian is chosen -% -% Output: -% scalo scalogram (squared magnitude of WT) -% f frequency samples (geometrically sampled between FMAX -% and FMIN). -% T time samples -% a scale vector (geometrically sampled between -% 1 and FMAX/FMIN) -% wt coefficient of the wavelet transform. -% X-axis corresponds to time (uniformly sampled), Y-axis -% corresponds to frequency (or scale) samples -% (geometrically sampled between Fmin (resp. Fmax/Fmin) -% and Fmax (resp. 1) -% First row of WT corresponds to the lowest analyzed -% frequency. -% wavescaled when the analyzing wavelet is Morlet or Mexican hat, -% wavescaled = wave. For an aritrary band-pass analyzing -% function, wavescaled contains columnwise the (N) -% scaled version of it -% -% Example: S = altes(256,0.1,0.45,10000) ; -% [scalo,f,T,a,wt,wavescaled] = contwtgn(S,0.01,0.5,128,8) ; -% -% See also: -% - -% Permission is granted for use and non-profit distribution providing that this -% notice be clearly maintained. The right to distribute any portion for profit -% or as part of any commercial product is specifically reserved for the author. - - -% CHECK INPUT FORMATS - -[xr,xc] = size(x) ; -if xr ~= 1 & xc ~= 1 - error('1-D signals only') -elseif xc == 1 - x = x.' ; -end - -T = 1 : length(x) ; - -% DEFAULT VALUES - -nt = size(x,2) ; -if exist('wave') == 0 - wave = 2 ; -end - -if nargin == 1 - XTF = fft(fftshift(x)) ; - sp = (abs(XTF(1:nt/2))).^2 ; - f = linspace(0,0.5,nt/2+1) ; f = f(1:nt/2) ; - plot(f,sp) ; grid ; - xlabel('frequency'); - title('Analyzed Signal Spectrum') ; - fmin = input('lower frequency bound = ') ; - fmax = input('upper frequency bound = ') ; - N = input('Frequency samples = ') ; - fmin_s = num2str(fmin) ; fmax_s = num2str(fmax) ; - N_s = num2str(N) ; - disp(['frequency runs from ',fmin_s,' to ',fmax_s,' over ',N_s,' points']) ; -end -if nargin == 4 | nargin == 5 - if fmin >= fmax - error('fmax must be greater than fmin') ; - end -end - -f = logspace(log10(fmax),log10(fmin),N) ; -a = logspace(log10(1),log10(fmax/fmin),N) ; amax = max(a) ; - -for ptr = 1:N - ha = gaussn(f(ptr),wave) ; nha = (length(ha)-1)/2 ; - detail = conv(x,ha) ; - wt(ptr,:) = detail(nha+1:nha+nt) ; -end - -wavescaled = wave ; - - -%%%% pour etre compatible avec le format de donnees de TFTB %%%%% - -wt = flipud(wt) ; -a = flipud(a(:)) ; -f = flipud(f(:)) ; -scalo = (wt.*conj(wt)) ; - - -%%%% pour etre compatible avec le format de donnees de TFTB %%%%% - - -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA diff --git a/mfiles/contwtgnmir.m b/mfiles/contwtgnmir.m deleted file mode 100644 index 9f2931b..0000000 --- a/mfiles/contwtgnmir.m +++ /dev/null @@ -1,84 +0,0 @@ -function [scalo,f,T,a,wt,wavescaled] = contwtgnmir(x,fmin,fmax,N,wave); - -% [scalo,f,T,a,wt,wavescaled] = contwtgnmir(x,fmin,fmax,N,wave); -% Continuous wavelet transform of mirrored 1-D signals -% If x = [a b c e f] is the signal to analyzed, contwtmir.m runs contwt.m -% on the mirrored version XxX = [c b [a b c d e f] e d]. The number of -% mirrored samples depends on the analyzed scale and the wavelet length. -% See contwt.m for Inputs/Outputs arguments. -% -% USE AN ORDER "wave" DERIVATIVE OF THE GAUSSIAN - -% CHECK INPUT FORMATS - -[xr,xc] = size(x) ; -if xr ~= 1 & xc ~= 1 - error('1-D signals only') -elseif xc == 1 - x = x.' ; -end - -T = 1 : lenfth(x) ; - -% DEFAULT VALUES - -nt = size(x,2) ; -if exist('wave') == 0 - wave = 2 ; -end - -if nargin == 1 - XTF = fft(fftshift(x)) ; - sp = (abs(XTF(1:nt/2))).^2 ; - f = linspace(0,0.5,nt/2+1) ; f = f(1:nt/2) ; - plot(f,sp) ; grid ; - xlabel('frequency'); - title('Analyzed Signal Spectrum') ; - fmin = input('lower frequency bound = ') ; - fmax = input('upper frequency bound = ') ; - N = input('Frequency samples = ') ; - fmin_s = num2str(fmin) ; fmax_s = num2str(fmax) ; - N_s = num2str(N) ; - disp(['frequency runs from ',fmin_s,' to ',fmax_s,' over ',N_s,' points']) ; -end -if nargin == 4 | nargin == 5 - if fmin >= fmax - error('fmax must be greater than fmin') ; - end -end - -f = logspace(log10(fmax),log10(fmin),N) ; -a = logspace(log10(1),log10(fmax/fmin),N) ; amax = max(a) ; - -for ptr = 1:N - ha = gaussn(f(ptr),wave) ; nha = (length(ha)-1)/2 ; - nbmir = min(nt,nha) ; - x_mir = [x(nbmir:-1:2) x x(nt-1:-1:nt-nbmir+1)] ; - detail = conv(x_mir,ha) ; - wt(ptr,1:nt) = detail(nha + nbmir : nha + nbmir + nt -1 ) ; -end -wavescaled = wave ; - -%%%% pour etre compatible avec le format de donnees de TFTB %%%%% - -wt = flipud(wt) ; -a = flipud(a(:)) ; -f = flipud(f(:)) ; -scalo = (wt.*conj(wt)) ; - - -%%%% pour etre compatible avec le format de donnees de TFTB %%%%% - -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA diff --git a/mfiles/correlmx.m b/mfiles/correlmx.m deleted file mode 100644 index 976a1a2..0000000 --- a/mfiles/correlmx.m +++ /dev/null @@ -1,67 +0,0 @@ -function Rx=correlmx(x,p,Rxtype); -% Rx=correlmx(x,p,Rxtype) correlation matrix of a signal -% -% Rx : correlation matrix (p+1) x (p+1) -% x : analyzed signal -% p : last autocorrelation lag -% Rxtype : computation algorithm (default : 'fbhermitian') -% possible values : 'hermitian', 'fbhermitian', 'burg' or 'fbburg' -% -% example : -% -% N=100; sig=real(fmconst(N,0.1))+0.4*randn(N,1); -% Rx=correlmx(sig,2,'burg'); [v,d] = eig(Rx), acos(-0.5*v(2,1)/v(1,1))/(2*pi) -% Rx=correlmx(sig,2,'hermitian'); [v,d] = eig(Rx), acos(-0.5*v(2,1)/v(1,1))/(2*pi) - -% F. Auger, july 1998. - -if (nargin<2), - error('At least two parameters required'); -elseif (nargin==2), - Rxtype='fbhermitian'; -end; - -[L,xcol]=size(x); -if xcol>1, - error('x must be a column vector'); -elseif p>L, - error('L must be greater than p'); -elseif p<1, - error('p must be greater than 0'); -end; - -Rxtype=upper(Rxtype); -if strcmp(Rxtype,'HERMITIAN')|strcmp(Rxtype,'FBHERMITIAN'), - vector=x(p+1-(0:p)); Rx=conj(vector) * vector.'; - for t=p+2:L, - vector=x(t-(0:p)); Rx=Rx+conj(vector) * vector.'; - end; - Rx=Rx/(L-p); - -elseif strcmp(Rxtype,'BURG')|strcmp(Rxtype,'FBBURG'), - R0=sum(abs(x).^2)/L; % variance - Rpos=zeros(1,p); Rneg=zeros(1,p); - for n=1:p, - Rpos(n)=sum(x(n+1:L).*conj(x(1:L-n)))/(L-n); - Rneg(n)=sum(x(1:L-n).*conj(x(n+1:L)))/(L-n); - end; - Rx=toeplitz([R0 Rpos],[R0 Rneg]); -else error(['unknown algorithm name' Rxtype]); -end; - -if strcmp(Rxtype(1:2),'FB'), - Rx=0.5*(Rx+Rx'); -end; -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA diff --git a/mfiles/d2statio.m b/mfiles/d2statio.m deleted file mode 100644 index 494d361..0000000 --- a/mfiles/d2statio.m +++ /dev/null @@ -1,42 +0,0 @@ -function [d,f]=d2statio(sig); -%D2STATIO Distance to stationarity -% [D,F]=D2STATIO(SIG) evaluates the distance of the signal -% to stationarity, using the pseudo Wigner-Ville distribution. -% -% SIG : signal to be analyzed (real or complex). -% D : vector giving the distance to stationarity for each frequency. -% F : vector of frequency bins -% -% Example : -% sig=noisecg(128); [d,f]=d2statio(sig); plot(f,d); -% xlabel('Frequency'); ylabel('Distance'); -% -% sig=fmconst(128); [d,f]=d2statio(sig); plot(f,d); -% xlabel('Frequency'); ylabel('Distance'); -% - -% O. Lemoine - May 1996. -% Copyright (c) by CNRS France, 1996. -% send bugs to f.auger@ieee.org - -N=length(sig); - -[tfr,t,f]=tfrspwv(sig); - -d2=((tfr-mean(tfr')'*ones(1,N))/norm(sig)).^2; - -d=mean(d2')'; - -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA diff --git a/mfiles/fmt.m b/mfiles/fmt.m deleted file mode 100644 index 331ec4d..0000000 --- a/mfiles/fmt.m +++ /dev/null @@ -1,138 +0,0 @@ -function [mellin,beta]=fmt(X,fmin,fmax,N); -%FMT Fast Fourier Mellin Transform. -% [MELLIN,BETA]=FMT(X,FMIN,FMAX,N) computes the Fast Mellin -% Transform of signal X. -% -% X : signal in time (Nx=length(X)). -% FMIN,FMAX : respectively lower and upper frequency bounds of -% the analyzed signal. These parameters fix the equivalent -% frequency bandwidth (expressed in Hz). When unspecified, you -% have to enter them at the command line from the plot of the -% spectrum. FMIN and FMAX must be >0 and <=0.5. -% N : number of analyzed voices. N must be even -% (default : automatically determined). -% MELLIN : the N-points Mellin transform of signal S. -% BETA : the N-points Mellin variable. -% -% Examples : -% sig=altes(128,0.05,0.45); -% [MELLIN,BETA]=fmt(sig,0.05,0.5,128); -% plot(BETA,real(MELLIN)); -% -% See also IFMT, FFT, IFFT. - -% P. Goncalves 9-95 - O. Lemoine, June 1996. -% Copyright (c) 1995 Rice University - CNRS (France) 1996. -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least one parameter required'); -end; - -[xrow,xcol] = size(X); - -if (nargin==2), - disp('FMIN will not be taken into account. Determine it with FMAX'); - disp(' from the following plot of the spectrum.'); -elseif nargin==3, - N=[]; -elseif (nargin==4 & rem(N,2)~=0), - error('N must be even'); -end; -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -Mt = length(X); -Z = hilbert(real(X)); -M = (Mt+rem(Mt,2))/2; - -if nargin<=2, % fmin,fmax,N unspecified - STF = fft(fftshift(X)); Nstf=length(STF); - sp = (abs(STF(1:Nstf/2))).^2; Maxsp=max(sp); - f = linspace(0,0.5,Nstf/2+1) ; f = f(1:Nstf/2); - plot(f,sp) ; grid; - xlabel('Normalized frequency'); - title('Analyzed signal energy spectrum'); - indmin=min(find(sp>Maxsp/1000)); - indmax=max(find(sp>Maxsp/1000)); - fmindflt=max([0.001 0.05*fix(f(indmin)/0.05)]); - fmaxdflt=0.05*ceil(f(indmax)/0.05); - txtmin=['Lower frequency bound [',num2str(fmindflt),'] : ']; - txtmax=['Upper frequency bound [',num2str(fmaxdflt),'] : ']; - fmin = input(txtmin); fmax = input(txtmax); - if fmin==[], fmin=fmindflt; end - if fmax==[], fmax=fmaxdflt; end -end - -if fmin >= fmax - error('FMAX must be greater or equal to FMIN'); -elseif fmin<=0.0 | fmin>0.5, - error('FMIN must be > 0 and <= 0.5'); -elseif fmax<=0.0 | fmax>0.5, - error('FMAX must be > 0 and <= 0.5'); -end - -B = fmax-fmin; % Bandwidth of the signal X -R = B/((fmin+fmax)/2); % Relative bandwidth of X - -Nq= ceil((B*Mt*(1+2/R)*log((1+R/2)/(1-R/2)))/2); -Nmin = Nq-rem(Nq,2); -Ndflt = 2^nextpow2(Nmin); -if nargin<=2, - Ntxt=['Number of frequency samples (>=',num2str(Nmin),') [',num2str(Ndflt),'] : ']; - N = input(Ntxt); -end - -if N~=[], - if (N= ',num2str(Nmin)]; - disp(dispstr); - end -else - N=Ndflt; -end - - -% Geometric sampling of the analyzed spectrum -No2 = N/2; -k = (1:No2); -q = (fmax/fmin)^(1/(No2-1)); -t = (1:Mt)-M-1; -geo_f = fmin*(exp((k-1).*log(q))); -tfmatx = zeros(Mt,N); -tfmatx = exp(-2*j*pi*t'*geo_f); -ZS = Z.'*tfmatx; -ZS(No2+1:N) = zeros(1,N-No2); - - -% Mellin transform computation of the analyzed signal -p = 0:(N-1); -L = log(fmin)/log(q); -mellin = N*log(q)*fftshift(ifft(ZS)).*exp(j*2*pi*L*(p/N-1/2)); -beta = (p/N-1/2)./log(q); - - -% Normalization -SP = fft(hilbert(real(X))); -indmin = 1+round(fmin*(xrow-2)); -indmax = 1+round(fmax*(xrow-2)); -SPana = SP(indmin:indmax); -nu = (indmin:indmax)'/N; -SPp = SPana./nu; -Normsig = sqrt(SPp'*SPana); - -mellin = mellin*Normsig/norm(mellin); diff --git a/mfiles/gaussn.m b/mfiles/gaussn.m deleted file mode 100644 index 03a84fa..0000000 --- a/mfiles/gaussn.m +++ /dev/null @@ -1,65 +0,0 @@ -function g = gaussn(f0,n) - -% function gn = gaussn(f0,n) : generates the order n derivative of the -% gaussian window, centered at frequency f0 -% The wavelet gn is real, but it is its analytic form that is -% synthetized first. The real part is then normalized by the analytical -% value of its energy (so, L2 normalization!) - -% alpha : std of the initial Gaussian -% g_0(t) = sqrt(2) (alpha/2*pi)^(1/4) exp(-alpha t^2) - -alpha = 2*pi^2*f0^2 / n ; - -% N = number of points, depends on "n". The support of the wavelet -% at tolerance 10^(-3) is determined experimentaly and modeled by -% a 4th order polynome (with reference frequency f0 is 0.1) - -PolyCoef = [-0.00000001420054 0.00001199042535 -0.00403466080616 ... - 0.92639865727446 20.64720217778273] ; -N = ceil((0.1/f0).*polyval(PolyCoef,n)) ; - -f = linspace(0,1,2*N+1) ; f = f(1:2*N) ; - -% Synthesis of the "analytic" wavelet Gn in the frequency domain - -G = (2*pi/alpha)^(1/4) * (2*i*pi*f).^n .* exp(-(pi^2*f.^2)./alpha) ; - -% Calculus of the wavelet energy (theoretical expression, cf. Gratshteyn, -% sec. 3.461, form. 2) - -p = (2*pi^2)/alpha ; -E = (2*pi/alpha)^(1/2) * (2*pi).^(2*n) * ... - (fact2(2*n-1))/(2*(2*p)^n)*sqrt(pi/p) ; -if isinf(E) | isnan(E) | E == 0 - E = integ(abs(G).^2,f) ; -end - -% L2 Normalized wavelet in time domain - -g = real ( fftshift(ifft(G))./(sqrt(E/2)) ) ; -g = g(2:2*N) ; -t = -(N-1):(N-1) ; - -% subplot(211) ; -% plot(f,abs(G)) ; title ('Wavelet Spectrum') ; -% xlabel('Frequency') ; grid -% subplot(212) -% plot(t,g) -% title('Wavelet in time') -% xlabel('Time') ; grid - - -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA diff --git a/mfiles/griffitc.m b/mfiles/griffitc.m deleted file mode 100644 index d63af0a..0000000 --- a/mfiles/griffitc.m +++ /dev/null @@ -1,46 +0,0 @@ -function [sig,iflaws]=griffitc(N,SNR); -%GRIFFITC Test signal example C of Griffiths' paper. -% [SIG,IFLAWS]=GRIFFITC(N,SNR) generates the test signal of -% the example C of the paper of Griffiths. -% -% N : signal length (default: 200) -% SNR : signal to noise ratio (default: 25 dB) -% SIG : output signal -% IFLAWS : instantaneous frequency laws of the 3 components -% -% Example : -% [sig,iflaws]=griffitc; plotifl(1:200,iflaws); grid; - -% F. Auger, July 1995. -% Copyright (c) 1996 by CNRS (France). -% Ref: L.J. Griffiths, "Rapid measurement of digital instantaneous -% frequency", IEEE Trans on ASSP, Vol 23, No 2, pp 207-222, 1975. -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin==0), - N=200; SNR=25; -elseif (nargin==1), - SNR=25; -end; -[sig1,iflaw1]=fmsin(N,0.25-0.08,0.25+0.08,192.6,50,0.285,+1); -[sig2,iflaw2]=fmsin(N,0.28-0.03,0.28+0.03,110.6,50,0.294,-1); -[sig3,iflaw3]=fmsin(N,0.40-0.02,0.40+0.02,149.6,50,0.417,-1); -noise=hilbert(randn(N,1)); -sig=sigmerge(sig1+sig2+sig3,noise,SNR); -if (nargout==2), - iflaws=[iflaw1 iflaw2 iflaw3]; -end; - diff --git a/mfiles/holder.m b/mfiles/holder.m deleted file mode 100644 index 7800f09..0000000 --- a/mfiles/holder.m +++ /dev/null @@ -1,99 +0,0 @@ -function h=holder(tfr,f,n1,n2,t,pl) -%HOLDER Estimate the Holder exponent through an affine TFR. -% H=HOLDER(TFR,F,N1,N2,T) estimates the Holder exponent of a -% function through an affine time-frequency representation of it, -% and plots the frequency marginal and the regression line. -% -% TFR : affine time-frequency representation. -% F : frequency values of the spectral analysis. -% N1 : indice of the minimum frequency for the linear regression. -% (default : 1). -% N2 : indice of the maximum frequency for the linear regression. -% (default : length(F)). -% T : time vector. If T is omitted, the function returns the -% global estimate of the Holder exponent. Otherwise, it -% returns the local estimates H(T) at the instants specified -% in T. -% H : output value (if T omitted) or vector (otherwise) containing -% the Holder estimate(s). -% -% Example : -% S=altes(128); [TFR,T,F]=tfrscalo(S,1:128,8); -% H=holder(TFR,F,1,length(F)); - -% P. Goncalves, October 1995 -% Copyright (c) 1995 Rice University -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -clf; -tfr = abs(tfr) ; - -if nargin<2, - error('There must be at least 2 input parameters'); -elseif nargin==2, - n1=1; n2=length(f); -elseif nargin==3, - n2=length(f); -end -if nargin<=5, - pl=1; -else - pl=0; -end - -if nargin <= 4 - fmarg = mean(tfr.').' ; - if pl, plot(log(f),log(fmarg)) , hold on; end; - p = polyfit(log(f(n1:n2)),log(fmarg(n1:n2)),1) ; - drt = p(1)*log(f)+p(2) ; - if pl, - plot(log(f),drt,'--g') ; - legend('Frequency marginal','Regression line'); - plot(log(f([n1 n2])),log(fmarg([n1 n2])),'+r') ; - xticklabels = round(logspace(log10(min(f)),log10(max(f)),4)*1000)./1000 ; - set(gca,'XTick',log(xticklabels)) ; - xlabel('frequency (logarithmically spaced)') ; - hold off ; - end - h = (-p(1)-1)/2 ; -else - if length(t) == 1 - if pl, plot(log(f),log(tfr(:,t))) ; hold on; end - p = polyfit(log(f(n1:n2)),log(tfr(n1:n2,t)),1) ; - drt = p(1)*log(f)+p(2) ; - if pl, - plot(log(f),drt,'--g') ; - legend('Frequency marginal','Regression line'); - plot(log(f([n1 n2])),log(tfr([n1 n2],t)),'+r') ; - xticklabels = round(logspace(log10(min(f)),log10(max(f)),4)*1000)./1000 ; - set(gca,'XTick',log(xticklabels)) ; - xlabel('frequency (logarithmically spaced)') ; - hold off ; - end - h = (-p(1)-1)/2 ; - elseif length(t) > 1 - [yt,xt] = size(t) ; if yt>xt, t = t.' ; end ; - j = 1 ; - for k = t, - p = polyfit(log(f(n1:n2)),log(tfr(n1:n2,k)),1) ; - h(j) = (-p(1)-1)/2 ; j = j + 1 ; - end - if pl, - plot(t,h) ;grid - title('Holder estimates at time instants T'); - end - end, -end diff --git a/mfiles/if2phase.m b/mfiles/if2phase.m deleted file mode 100644 index f947ea8..0000000 --- a/mfiles/if2phase.m +++ /dev/null @@ -1,26 +0,0 @@ -function phi=if2phase(iflaw) -%IF2PHASE Generate the phase from the instantaneous frequency - -% O. Lemoine - February 1996. -% experimental - -N=length(iflaw); -phi=zeros(N,1); - -%phi=2*pi*cumsum(iflaw); -for k=1:N, - t=1:k; - phi(k)=2*pi*integ(iflaw(t),t); -end% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA diff --git a/mfiles/ifestar2.m b/mfiles/ifestar2.m deleted file mode 100644 index 1254765..0000000 --- a/mfiles/ifestar2.m +++ /dev/null @@ -1,91 +0,0 @@ -function [fnorm,t,rejratio]=ifestar2(x,t); -%IFESTAR2 Instantaneous frequency estimation using AR2 modelisation. -% [FNORM,T2,RATIO]=IFESTAR2(X,T) computes an estimate of the -% instantaneous frequency of the real signal X at time -% instant(s) T. The result FNORM lies between 0.0 and 0.5. This -% estimate is based only on the 4 last signal points, and has -% therefore an approximate delay of 2.5 points. -% -% X : real signal to be analyzed. -% T : Time instants (must be greater than 4) -% (default : 4:length(X)). -% FNORM : Output (normalized) instantaneous frequency. -% T2 : Time instants coresponding to FNORM. Since the -% algorithm can not always give a value, T2 is -% different of T. -% RATIO : proportion of instants where the algorithm yields -% an estimation -% -% Examples : -% [x,if]=fmlin(50,0.05,0.3,5); x=real(x); [if2,t]=ifestar2(x); -% plot(t,if(t),t,if2); -% -% N=1100; [deter,if]=fmconst(N,0.05); deter=real(deter); -% noise=randn(N,1); NbSNR=101; SNR=linspace(0,100,NbSNR); -% for iSNR=1:NbSNR, -% sig=sigmerge(deter,noise,SNR(iSNR)); -% [if2,t,ratio(iSNR)]=ifestar2(sig); -% EQM(iSNR)=norm(if(t)-if2)^2 / length(t) ; -% end; -% subplot(211); plot(SNR,-10.0 * log10(EQM)); grid; -% xlabel('SNR'); ylabel('-10 log10(EQM)'); -% subplot(212); plot(SNR,ratio); grid; -% xlabel('SNR'); ylabel('ratio'); -% -% See also INSTFREQ, KAYTTH, SGRPDLAY. - -% F. Auger, April 1996. -% This estimator is the causal version of the estimator called -% "4 points Prony estimator" in the article "Instantaneous -% frequency estimation using linear prediction with comparisons -% to the dESAs", IEEE Signal Processing Letters, Vo 3, No 2, p -% 54-56, February 1996. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least one parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol~=1), - error('X must have only one column'); -end -x = real(x); - -if (nargin == 1), - t=4:xrow; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('T must have only one row'); -elseif min(t)<4, - error('The smallest value of T must be greater than 4'); -end; - - -Kappa = x(t-1) .* x(t-2) - x(t ) .* x(t-3) ; -psi1 = x(t-1) .* x(t-1) - x(t ) .* x(t-2) ; -psi2 = x(t-2) .* x(t-2) - x(t-1) .* x(t-3) ; -den = psi1 .* psi2 ; -indices = find(den>0); -arg=0.5*Kappa(indices)./sqrt(den(indices)); -indarg=find(abs(arg)>1); -arg(indarg)=sign(arg(indarg)); -fnorm = acos(arg)/(2.0*pi); -rejratio = length(indices)/length(t); -t = t(indices); - diff --git a/mfiles/ifmt.m b/mfiles/ifmt.m deleted file mode 100644 index 015e879..0000000 --- a/mfiles/ifmt.m +++ /dev/null @@ -1,86 +0,0 @@ -function x=ifmt(mellin,beta,M); -%IFMT Inverse fast Mellin transform. -% X=IFMT(MELLIN,BETA,M) computes the inverse fast Mellin -% transform of MELLIN. -% WARNING : the inverse of the Mellin transform is correct only -% if the Mellin transform has been computed from FMIN to 0.5 Hz, -% and if the original signal is analytic. -% -% MELLIN : Mellin transform to be inverted. Mellin must have been -% obtained from FMT with frequency running from FMIN to 0.5 Hz. -% BETA : Mellin variable issued from FMT. -% M : number of points of the inverse Mellin transform. -% (default : length(MELLIN)). -% X : inverse Mellin transform with M points in time. -% -% Example : -% sig=atoms(128,[64,0.25,32,1]); -% [MELLIN,BETA]=fmt(sig,0.05,0.5,256); clf; -% X=ifmt(MELLIN,BETA,128); plot(real(X)); hold; -% plot(real(sig),'g'); hold; -% -% See also : fmt, fft, ifft. -% - -% P. Goncalves 9-95 - O. Lemoine, June 1996. -% Copyright (c) 1995 Rice University - CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -N=length(mellin); - -if nargin<2, - error('At least 2 input parameters required'); -elseif nargin==2, - M=N; -end - -No2 = (N+rem(N,2))/2; -q = exp(1/(N*(beta(2)-beta(1)))); -fmin = 0.5/(q^(No2-1)); - - -% Inverse Mellin transform computation -p = 0:(N-1); -L = log(fmin)/log(q); -S = fft(fftshift(mellin.*exp(-j*2*pi*L*(p/N-1/2))/(N*log(q)))); -S = S(1:No2); - - -% Inverse Fourier transform -k = (1:No2); -x = zeros(M,1); -t = (1:M)-(M+rem(M,2))/2-1; -geo_f = fmin*(exp((k-1).*log(q))) ; -itfmatx = zeros(M,No2); -itfmatx = exp(2*i*t'*geo_f*pi); -for k=1:M, - x(k) = real(integ(itfmatx(k,:).*S,geo_f)); -end; - -x = hilbert(x); - -% Normalization -SP = fft(x); fmax=0.5; -indmin = 1+round(fmin*(M-2)); -indmax = 1+round(fmax*(M-2)); -SPana = SP(indmin:indmax); -nu = (indmin:indmax)'/M; -SPp = SPana./nu; -Esm = SPp'*SPana; -x = x*norm(mellin)/sqrt(Esm); - - - diff --git a/mfiles/imextrac.m b/mfiles/imextrac.m deleted file mode 100644 index 5d5cd87..0000000 --- a/mfiles/imextrac.m +++ /dev/null @@ -1,81 +0,0 @@ -function [Image2,NbDots]=imextrac(Image,trace); -% imextrac(Image) extract and isolate dots in a binary image -% -% example: -% Image=[1 0 0 0 0 0 1 0 0 0 1 0 ;... -% 1 1 0 0 0 1 1 0 0 0 1 1 ;... -% 1 0 0 0 1 0 1 0 0 0 0 1 ;... -% 0 0 0 0 1 1 1 0 0 0 0 0 ;... -% 0 0 0 0 0 0 0 0 0 0 0 0 ]; -% image2=imextrac(Image) -% -% See also tfrsurf. - -% F. Auger, oct 1999 -% Copyright (c) CNRS - France 1999. -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -if nargin==1, trace=0; end; - -Image2=Image; -[Nbrow,Nbcol]=size(Image2); - -NbDots=1 ; -if trace==1, fprintf('extracting dots in the image.\n'); end; - -if (Image2(1,1)==1), - NbDots=NbDots+1; Image2(1,1)=NbDots; -end; - -for i=2:Nbcol, - if (Image2(1,i)==1), - if (Image2(1,i-1)>1), - Image2(1,i)=Image2(1,i-1); - else - NbDots=NbDots+1; Image2(1,i)=NbDots; - end; - end; -end; -if trace==1, disprog(1,Nbrow,10); end; - -for j=1:Nbrow, - if (Image2(j,1)==1), - if (Image2(j-1,1)>1), - Image2(j,1)=Image2(j-1,1); - else - NbDots=NbDots+1; Image2(j,1)=NbDots; - end; - end; - - for i=2:Nbcol, - if (Image2(j,i)==1), - if (Image2(j-1,i)==0)&(Image2(j,i-1)==0), - NbDots=NbDots+1; Image2(j,i)=NbDots; - elseif (Image2(j-1,i)==0)&(Image2(j,i-1)>1), - Image2(j,i)=Image2(j,i-1); - elseif (Image2(j-1,i)>1)&(Image2(j,i-1)==0), - Image2(j,i)=Image2(j-1,i); - elseif (Image2(j-1,i)>1)&(Image2(j,i-1)>1), - MinDot=min(Image2(j-1,i),Image2(j,i-1)); - MaxDot=max(Image2(j-1,i),Image2(j,i-1)); - Indices=find(Image2==MaxDot); Image2(Indices)=MinDot; - Image2(j,i)=MinDot; - else error('should never happen'); - end; - end; - end; - if trace==1, disprog(j,Nbrow,10); end; -end; - diff --git a/mfiles/izak.m b/mfiles/izak.m deleted file mode 100644 index 507253c..0000000 --- a/mfiles/izak.m +++ /dev/null @@ -1,42 +0,0 @@ -function sig=izak(dzt) -%IZAK Inverse Zak transform. -% SIG=IZAK(DZT) computes the inverse Zak transform of matrix DZT. -% -% DZT : (N,M) matrix of Zak samples. -% SIG : Output signal (M*N,1) containing the inverse Zak transform. -% -% Example : -% sig=fmlin(256); DZT=zak(sig); sigr=izak(DZT); -% plot(real(sigr)); hold; plot(real(sig)); hold; -% -% See also ZAK, TFRGABOR. - -% O. Lemoine - February 1996 -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -[N,M]=size(dzt); - -if nargin<1, - error('The number of parameters must be one'); -end - -sig=zeros(N*M,1); - -for m=1:M, - sig(m+(0:N-1)*M)=sqrt(N)*ifft(dzt(:,m)); -end - diff --git a/mfiles/kaytth.m b/mfiles/kaytth.m deleted file mode 100644 index 45b222d..0000000 --- a/mfiles/kaytth.m +++ /dev/null @@ -1,29 +0,0 @@ -function H=kaytth(length); -%KAYTTH Kay-Tretter filter computation. -% H=KAYTTH(length); Kay-Tretter filter computation. -% -% See also INSTFREQ - -% F. Auger, March 1994. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -pp1=length*(length+1); -den=2.0*length*(length+1)*(2.0*length+1.0)/3.0; -i=1:length; H=pp1-i.*(i-1); - -H=H ./ den; - diff --git a/mfiles/margtfr.m b/mfiles/margtfr.m deleted file mode 100644 index 6e20f71..0000000 --- a/mfiles/margtfr.m +++ /dev/null @@ -1,50 +0,0 @@ -function [margt,margf,E]=margtfr(tfr,t,f) -%MARGTFR Marginals and energy of a time-frequency representation. -% [MARGT,MARGF,E]=MARGTFR(TFR,T,F) calculates the time and -% frequency marginals and the energy of a time-frequency -% representation. -% -% TFR : time-frequency representation (M,N) -% T : vector containing the time samples in sec. -% (default : (1:N)) -% F : vector containing the frequency samples in Hz, not -% necessary uniformly sampled. (default : (1:M)) -% MARGT : time marginal -% MARGF : frequency marginal -% E : energy of TFR -% -% Example : -% S=altes(128,0.05,0.45); TFR=tfrscalo(S,1:128,8,'auto'); -% [MARGT,MARGF,E] = margtfr(TFR); -% subplot(211); plot(T,MARGT); subplot(212); plot(F,MARGF); -% -% See also MOMTTFR. - -% P. Goncalves, October 95 -% Copyright (c) 1995 Rice University -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -[M,N] = size(tfr) ; -if nargin == 1, - t = (1:N); - f = (1:M)'; -elseif nargin == 2, - f = (1:M)'; -end - -E = real(integ2d(tfr,t,f)/M) ; -margt = real(integ(tfr.',f')/M) ; -margf = real(integ(tfr,t)/N) ; diff --git a/mfiles/midpoint.m b/mfiles/midpoint.m deleted file mode 100644 index dc878f7..0000000 --- a/mfiles/midpoint.m +++ /dev/null @@ -1,87 +0,0 @@ -function [ti,fi]=midpoint(t1,f1,t2,f2,k) -%MIDPOINT Mid-point construction used in the interference diagram. -% [TI,FI]=MIDPOINT(T1,F1,T2,F2,K) gives the coordinates in the -% time-frequency plane of the interference-term corresponding to -% the points (T1,F1) and (T2,F2), for a distribution in the -% affine class perfectly localized on power-law group-delays of -% the form : tx(nu)=t0+c nu^(K-1). -% -% T1 : time-coordinate of the first point -% F1 : frequency-coordinate of the first point (>0) -% T2 : time-coordinate of the second point -% F2 : frequency-coordinate of the second point (>0) -% K : power of the group-delay law -% K = 2 : Wigner-Ville -% K = 1/2 : D-Flandrin -% K = 0 : Bertrand (unitary) -% K = -1 : Unterberger (active) -% K = inf : Margenau-Hill-Rihaczek -% TI : time-coordinate of the interference term -% FI : frequency-coordinate of the interference term -% -% See also PLOTSID. - -% P. Flandrin, September 1995 - F. Auger, April 1996. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if f1<=0 | f2<=0, - error('F1 and F2 must be >0'); -end -[rt1,ct1]=size(t1); -[rt2,ct2]=size(t2); -[rf1,cf1]=size(f1); -[rf2,cf2]=size(f2); -if (rt1~=rt2|rt1~=rf1|rt1~=rf2) | (ct1~=ct2|ct1~=cf1|ct1~=cf2), - error('T1, T2, F1 and F2 must have the same size'); -end -if rt1>ct1, - error('T1 must be a row-vector'); -elseif rt2>ct2, - error('T2 must be a row-vector'); -elseif rf2>cf2, - error('F2 must be a row-vector'); -elseif rf1>cf1, - error('F1 must be a row-vector'); -end - -if (k==2), - fi=(f1+f2)/2; - ti=(t1+t2)/2; -elseif (k==inf), - ti=[t1;t2]; - fi=[f2;f1]; -else - I=find(abs(f1-f2)>sqrt(eps)); - if length(I)~=0, - if (k==1), - fi(I)=exp( (f1(I).*(log(f1(I))-1)-f2(I).*(log(f2(I))-1)) ./ ... - (f1(I)-f2(I))); - ti(I)=(t1(I).*f1(I)-t2(I).*f2(I)) ./ (f1(I)-f2(I)) - ... - (t1(I)-t2(I)) ./ (log(f1(I))-log(f2(I))); - elseif (k==0), - fi(I)=(f1(I)-f2(I))./(log(f1(I))-log(f2(I))); - ti(I)=(t1(I).*f1(I)-t2(I).*f2(I)) ./ (f1(I)-f2(I)) + ... - f1(I) .* f2(I) .* (t2(I)-t1(I)) .* ... - (log(f1(I))-log(f2(I))) ./ (f2(I)-f1(I)).^2; - else - t0(I)=(t1(I).*f2(I).^(k-1)-t2(I).*f1(I).^(k-1)) ./ ... - (f2(I).^(k-1)-f1(I).^(k-1)); - fi(I)=((f1(I).^k-f2(I).^k) ./ (f1(I)-f2(I))/k).^(1/(k-1)); - ti(I)=t0(I)+(t2(I)-t1(I)) ./ (f2(I).^(k-1)-f1(I).^(k-1)) .*fi(I).^(k-1); - end - end -end diff --git a/mfiles/momftfr.m b/mfiles/momftfr.m deleted file mode 100644 index b8f17e3..0000000 --- a/mfiles/momftfr.m +++ /dev/null @@ -1,56 +0,0 @@ -function [tm,D2]=momftfr(tfr,tmin,tmax,time); -%MOMFTFR Frequency moments of a time-frequency representation. -% [TM,D2]=MOMFTFR(TFR,TMIN,TMAX,TIME) computes the frequeny -% moments of a time-frequency representation. -% -% TFR : time-frequency representation ([Nrow,Ncol]size(TFR)). -% TMIN : smallest column element of TFR taken into account -% (default : 1) -% TMAX : highest column element of TFR taken into account -% (default : Ncol) -% TIME : true time instants (default : 1:Ncol) -% TM : averaged time (first order moment) -% D2 : squared time duration (second order moment) -% -% Example : -% sig=fmlin(128,0.1,0.4); -% [tfr,t,f]=tfrwv(sig); [tm,D2]=momftfr(tfr); -% subplot(211); plot(f,tm); subplot(212); plot(f,D2); -% -% See also MOMTTFR, MARGTFR. - -% F. Auger, August 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -[tfrrow,tfrcol]=size(tfr); -if (nargin==1), - tmin=1; tmax=tfrcol; time=tmin:tmax; -elseif (nargin==2), - tmax=tfrcol; time=tmin:tmax; -elseif (nargin==3), - time=tmin:tmax; -end; - -if (tmin>tmax)|(tmin<=0)|(tmax>tfrcol), - error('1<=TMIN<=TMAX<=Ncol'); -end; - -E = sum(tfr(:,tmin:tmax).'); -tm = (time * tfr(:,tmin:tmax).' ./E).'; -D2 = (time.^2 * tfr(:,tmin:tmax).' ./E).' - tm.^2; - - diff --git a/mfiles/momttfr.m b/mfiles/momttfr.m deleted file mode 100644 index 1c06e6f..0000000 --- a/mfiles/momttfr.m +++ /dev/null @@ -1,98 +0,0 @@ -function [fm,B2]=momttfr(tfr,method,fbmin,fbmax,freqs); -%MOMTTFR Time moments of a time-frequency representation. -% [FM,B2]=MOMTTFR(TFR,METHOD,FBMIN,FBMAX,FREQS) computes the time -% moments of a time-frequency representation. -% -% TFR : time-frequency representation ([Nrow,Ncol]=size(TFR)). -% METHOD: chosen representation (name of the corresponding M-file). -% FBMIN : smallest frequency bin (default : 1) -% FBMAX : highest frequency bin (default : Nrow) -% FREQS : true frequency of each frequency bin. FREQS must be of -% length FBMAX-FBMIN+1. -% (default : 0:step:(0.5-step) or -0.5:step:(0.5-step) -% depending on METHOD) -% FM : averaged frequency (first order moment) -% B2 : frequency band (second order moment) -% -% Example : -% sig=fmlin(128,0.1,0.4); tfr=tfrwv(sig); -% [fm,B2]=momttfr(tfr,'tfrwv'); -% subplot(211); plot(fm); subplot(212); plot(B2); -% freqs=linspace(0,63/128,64); tfr=tfrsp(sig); -% [fm,B2]=momttfr(tfr,'tfrsp',1,64,freqs); -% subplot(211); plot(fm); subplot(212); plot(B2); -% -% See also MOMFTFR, MARGTFR. - -% F. Auger, August 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -[tfrrow,tfrcol]=size(tfr); -if (nargin<=1), - error('At least two input arguments'); -elseif (nargin==2), - fbmin=1; fbmax=tfrrow; -elseif (nargin==3), - fbmax=tfrrow; -end -if (fbmin>fbmax)|(fbmin==0)|(fbmax>tfrrow), - error('1<=FBMIN<=FBMAX<=Nrow'); -elseif nargin==5, - if length(freqs)~=(fbmax-fbmin+1), - error('FREQS must have FBMAX-FBMIN+1 elements'); - end -end; - -method=upper(method); -if strcmp(method,'TFRWV' ) | strcmp(method,'TFRPWV' ) | ... - strcmp(method,'TFRSPWV' ) | strcmp(method,'TFRCW' ) | ... - strcmp(method,'TFRZAM' ) | strcmp(method,'TFRBJ' ) | ... - strcmp(method,'TFRBUD' ) | strcmp(method,'TFRGRD' ) | ... - strcmp(method,'TFRRSPWV') | strcmp(method,'TFRRPWV' ) | ... - strcmp(method,'TFRRIDB' ) | strcmp(method,'TFRRIDH' ) | ... - strcmp(method,'TFRRIDT' ) | strcmp(method,'TFRASPW' ) | ... - strcmp(method,'TFRDFLA' ) | strcmp(method,'TFRSPBK' ) | ... - strcmp(method,'TFRUNTAC') | strcmp(method,'TFRUNTPA') | ... - strcmp(method,'TFRBERT' ) | strcmp(method,'TFRSCALO') | ... - strcmp(method,'TYPE2' ), - typertf='TYPE2'; -elseif strcmp(method,'TFRPMH' )| strcmp(method,'TFRRPMH' )| ... - strcmp(method,'TFRSP' )| strcmp(method,'TFRRSP' )| ... - strcmp(method,'TFRPPAGE')| strcmp(method,'TFRRPPAG')| ... - strcmp(method,'TFRMHS' )| strcmp(method,'TFRRGAB' )| ... - strcmp(method,'TFRMH' )| strcmp(method,'TFRMMCE' )| ... - strcmp(method,'TFRMSC' )| strcmp(method,'TFRRMSC' )| ... - strcmp(method,'TFRPAGE' )| strcmp(method,'TFRGABOR')| ... - strcmp(method,'TFRRI' )| strcmp(method,'TYPE1' ), - typertf='TYPE1'; -else - error('Unknown representation.'); -end; - -if nargin<=4, - if strcmp(typertf,'TYPE1'), - freqs=rem((fbmin-1:fbmax-1)/tfrrow + 0.5, 1.0) - 0.5; - elseif strcmp(typertf,'TYPE2'), - freqs=0.5*(fbmin-1:fbmax-1)/tfrrow; - end; -end - -E = sum(tfr(fbmin:fbmax,:)); -fm = (freqs * tfr(fbmin:fbmax,:) ./E).'; -B2 = (freqs.^2 * tfr(fbmin:fbmax,:) ./E).' - fm.^2; - - diff --git a/mfiles/parafrep.m b/mfiles/parafrep.m deleted file mode 100644 index 1a3b3d9..0000000 --- a/mfiles/parafrep.m +++ /dev/null @@ -1,92 +0,0 @@ -function [spec,freqs,A]=parafrep(Rx,N,method,q); -% [spec,freqs]=parafrep(Rx,N,method,q) parametric frequency representation -% of a signal. -% -% Rx : correlation matrix of size (p+1)x(p+1) -% N : number of frequency bins between 0 and 0.5 -% method : can be either 'ar', 'periodogram', 'capon', 'capnorm', 'lagunas', -% or 'genlag'. -% q : parameter for the generalized Lagunas method. -% -% noise=rand(1000,1); signal=filter([1 0 0],[1 1 1],noise); -% figure(1);parafrep(correlmx(signal,2,'hermitian'),128,'AR');title('AR (2)'); -% figure(2);parafrep(correlmx(signal,4,'hermitian'),128,'Capon');title('Capon (4)'); -% figure(3);parafrep(correlmx(signal,2,'hermitian'),128,'lagunas');title('Lagunas (2)'); -% figure(4);parafrep(correlmx(signal,40,'hermitian'),128,'periodogram');title('periodogram (40)'); - -% F. Auger, july 1998, april 99. - -if nargin<1, - error('At least one parameter required'); -elseif nargin==1, - N=128; method='capon'; -elseif nargin==2, - method='capon'; -end; - -[Rxrow,Rxcol]=size(Rx); -if (Rxrow ~= Rxcol), - error('Rx must be a square matrix'); -end; - -p=Rxrow-1; -freqs=linspace(0,0.5,N); -spec=zeros(N,1); - -method=upper(method); -if strcmp(method,'AR'), - Un=ones(Rxrow,1); Rxm1Un= (Rx\Un); - P1=real(Un'*Rxm1Un); A=Rxm1Un/P1; - for ifreq=1:N, - Z=exp(2j*pi*freqs(ifreq)*(0:Rxrow-1)'); - spec(ifreq)=P1 ./ abs(Z' * A)^2 ; - end; -elseif strcmp(method,'PERIODOGRAM'), - for ifreq=1:N, - Z=exp(2j*pi*freqs(ifreq)*(0:Rxrow-1)'); - spec(ifreq)=real(Z' * Rx *Z)/(p+1)^2; - end; -elseif strcmp(method,'CAPON'), - for ifreq=1:N, - Z=exp(2j*pi*freqs(ifreq)*(0:Rxrow-1)'); - spec(ifreq)=1.0 / real(Z' * (Rx\Z)); - end; -elseif strcmp(method,'CAPNORM'), - for ifreq=1:N, - Z=exp(2j*pi*freqs(ifreq)*(0:Rxrow-1)'); - spec(ifreq)=(p+1) / real(Z' * (Rx\Z)); - end; -elseif strcmp(method,'LAGUNAS'), - for ifreq=1:N, - Z=exp(2j*pi*freqs(ifreq)*(0:Rxrow-1)'); - Rxm1Z=Rx\Z; spec(ifreq)=real(Z' * Rxm1Z)/real(Z' * (Rx\Rxm1Z)); - end; -elseif strcmp(method,'GENLAG'), - for ifreq=1:N, - Z=exp(2j*pi*freqs(ifreq)*(0:Rxrow-1)'); - Rxqm1Z=(Rx)^q \Z; spec(ifreq)=real(Z' * Rx * Rxqm1Z)/real(Z' * (Rx\Rxqm1Z)); - end; -else - error('unknown frequency representation'); -end; - -if (nargout==0), - figure(gcf); plot(freqs,10.0*log10(spec)); grid; - xlabel('normalized frequency'); - ylabel('DSP (dB)'); -end; - - -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA diff --git a/mfiles/plotsid.m b/mfiles/plotsid.m deleted file mode 100644 index 1c8ee9a..0000000 --- a/mfiles/plotsid.m +++ /dev/null @@ -1,107 +0,0 @@ -function plotsid(t,iflaws,k); -%PLOTSID Schematic interference diagram of FM signals. -% PLOTSID(T,IFLAWS,K) plots the schematic interference diagram of -% (analytic) FM signals. -% -% T : time instants, -% IFLAWS : matrix of instantaneous frequencies, -% with as may columns as signal components. -% K : distribution (default : 2): -% K = 2 : Wigner-Ville -% K = 1/2 : D-Flandrin -% K = 0 : Bertrand (unitary) -% K = -1 : Unterberger (active) -% K = inf : Margenhau-Hill-Rihaczek -% -% Example : -% Nt=90; [y,iflaw]=fmlin(Nt,0.05,0.25); -% [y2,iflaw2]=fmconst(50,0.4); -% iflaw(:,2)=[NaN*ones(10,1);iflaw2;NaN*ones(Nt-60,1)]; -% plotsid(1:Nt,iflaw,0); -% -% See also PLOTIFL, MIDPOINT. - -% P. Flandrin, September 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin==1), - error('at least 2 parameters : t and iflaws'); -elseif (nargin==2), - k=2; % Wigner-Ville -end; - -indices=find(1-isnan(iflaws)); -if (min(iflaws(indices))<0)|(max(iflaws(indices))>0.5), - error ('each element of IFLAWS must be between 0 and 0.5'); -end; - -[iflawrow,iflawcol]=size(iflaws); -tcol=length(t); -clf; figure(gcf); plotifl(t,iflaws); -hold on; -col=['y','m','c','r']; - -% auto-terms -for j=1:iflawcol, - indices=find(1-isnan(iflaws(:,j))); - Nbpoints=length(indices); - for i=1:Nbpoints-1, - ta= t(indices(i)) *ones(1,Nbpoints-i); - fa= iflaws(indices(i),j)*ones(1,Nbpoints-i); - tb= t(indices(i+1:Nbpoints)); - fb=iflaws(indices(i+1:Nbpoints),j)'; - [ti,fi]=midscomp(ta,fa,tb,fb,k); - plot(ti,fi,['.',num2str(col(rem(j-1,4)+1))]); - end; -end; - -% cross-terms -for j1=1:iflawcol, - indices1=find(1-isnan(iflaws(:,j1))); - Nbpoints1=length(indices1); - for j2=j1+1:iflawcol, - indices2=find(1-isnan(iflaws(:,j2))); - Nbpoints2=length(indices2); - for i=1:Nbpoints1, - ta= t(indices1(i)) *ones(1,Nbpoints2); - fa= iflaws(indices1(i),j1)*ones(1,Nbpoints2); - tb= t(indices2); - fb=iflaws(indices2,j2)'; - [ti,fi]=midscomp(ta,fa,tb,fb,k); - plot(ti,fi,'.g') - end; - end; -end; - -hold off -axis([t(1) t(tcol) 0 0.5]); -grid -if k==2, - dist=' of the Wigner-Ville distribution'; -elseif k==1/2, - dist=' of the D-Flandrin distribution'; -elseif k==0, - dist=' of the (unitary) Bertrand distribution'; -elseif k==-1, - dist=' of the (active) Unterberger distribution'; -elseif k>1/sqrt(eps), - dist=' of the Margenhau-Hill-Rihaczek distribution'; -else - dist=''; -end - -title(['Interference diagram',dist,' (k = ',num2str(k),')']); diff --git a/mfiles/scale.m b/mfiles/scale.m deleted file mode 100644 index 91dd30d..0000000 --- a/mfiles/scale.m +++ /dev/null @@ -1,151 +0,0 @@ -function S=scale(X,a,fmin,fmax,N,trace); -%SCALE Scale a signal using the Mellin transform. -% S=SCALE(X,A,FMIN,FMAX,N,TRACE) computes the A-scaled version -% of signal X : A^(-1/2) X(T/A) using its Mellin transform. -% -% X : signal in time to be scaled (Nx=length(X)). -% A : scale factor. A < 1 corresponds to a compression in the time -% domain. A can be a vector. (default : 2) -% FMIN,FMAX : respectively lower and upper frequency bounds of -% the analyzed signal. These parameters fix the equivalent -% frequency bandwidth (expressed in Hz). When unspecified, you -% have to enter them at the command line from the plot of the -% spectrum. FMIN and FMAX must be >0 and <=0.5. -% N : number of analyzed voices (default : automatically determined). -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% S : the A-scaled version of signal X. Length of S can be larger -% than length of X if A > 1. If A is a vector of length L, S is -% a matrix with L columns. S has the same energy as X. -% -% Example : -% sig=klauder(128); S=scale(sig,2,.05,.45,128); -% subplot(211); plot(sig); subplot(212); plot(real(S(65:192))); - -% P. Goncalves, October 1995 - O. Lemoine, June 1996. -% Copyright (c) Rice University - CNRS (France) -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least one parameter required'); -end; - -[Mt,xcol] = size(X); - -if (nargin == 1), - a=2; -elseif (nargin==3), - disp('FMIN will not be taken into account. Determine it with FMAX'); - disp(' from the following plot of the spectrum.'); -elseif nargin==4, - N=[]; -end; - -if nargin <=5, - trace=0; -end -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end; - -Z = hilbert(real(X)); -T = Mt; -M = (Mt+rem(Mt,2))/2; - -if nargin<=3 % fmin,fmax,N unspecified - STF = fft(fftshift(Z)); Nstf=length(STF); - sp = (abs(STF(1:Nstf/2))).^2; Maxsp=max(sp); - f = linspace(0,0.5,Nstf/2+1) ; f = f(1:Nstf/2); - plot(f,sp) ; grid; - xlabel('Normalized frequency'); - title('Analyzed signal energy spectrum'); - indmin=min(find(sp>Maxsp/1000)); - indmax=max(find(sp>Maxsp/1000)); - fmindflt=max([0.01 0.05*fix(f(indmin)/0.05)]); - fmaxdflt=0.05*ceil(f(indmax)/0.05); - txtmin=['Lower frequency bound [',num2str(fmindflt),'] : ']; - txtmax=['Upper frequency bound [',num2str(fmaxdflt),'] : ']; - fmin = input(txtmin); fmax = input(txtmax); - if fmin==[], fmin=fmindflt; end - if fmax==[], fmax=fmaxdflt; end -end - -if fmin >= fmax - error('FMAX must be greater or equal to FMIN'); -elseif fmin<=0.0 | fmin>0.5, - error('FMIN must be > 0 and <= 0.5'); -elseif fmax<=0.0 | fmax>0.5, - error('FMAX must be > 0 and <= 0.5'); -end - -B = fmax-fmin ; -R = B/((fmin+fmax)/2) ; - -Nq= ceil((B*T*(1+2/R)*log((1+R/2)/(1-R/2)))/2); -Nmin = Nq-rem(Nq,2); -Ndflt = 2^nextpow2(Nmin); -if nargin<=3, - Ntxt=['Number of frequency samples (>=',num2str(Nmin),') [',num2str(Ndflt),'] : ']; - N = input(Ntxt); -end -if N~=[], - if (N ',num2str(Nmin)]; - disp(dispstr); - end -else - N=Ndflt; -end - - -% Geometric sampling of the analyzed spectrum -k = 1:N; -q = (fmax/fmin)^(1/(N-1)); -geo_f = fmin*(exp((k-1).*log(q))); -t = (1:Mt)-M-1; -tfmatx = zeros(Mt,N); -tfmatx = exp(-2*i*t'*geo_f*pi); -ZS = Z.'*tfmatx; -ZS(N+1:2*N) = zeros(1,N); - - -% Mellin transform computation of the analyzed signal -p = 0:(2*N-1); -MS = fftshift(ifft(ZS)); -beta = (p/N-1)./(2*log(q)); - - -% Inverse Mellin transform and inverse Fourier transform -Mmax = max(ceil(Mt/2*a)); -S = zeros(2*Mmax,length(a)); -ptr = 1; -for acurrent = a, - if trace, disprog(ptr,length(a),10); end - DMS = exp(-2*i*pi*beta*log(acurrent)).*MS; % Scaling in Mellin domain - DS = fft(fftshift(DMS)); % Inverse Mellin transform - Mcurrent = ceil(acurrent*Mt/2); - t = [-Mcurrent:Mcurrent-1]-1; - itfmatx = zeros(2*Mcurrent,N); - itfmatx = exp(2*i*t'*geo_f*pi); - dilate_sig = zeros(2*Mcurrent,1); - for kk=1:2*Mcurrent, - dilate_sig(kk) = integ(itfmatx(kk,:).*DS(1:N),geo_f) ; - end; - S(Mmax-Mcurrent+1:Mmax+Mcurrent,ptr) = dilate_sig; - ptr=ptr+1; -end - -S=S*norm(X)/norm(S); diff --git a/mfiles/tffilter.m b/mfiles/tffilter.m deleted file mode 100644 index d28a0fc..0000000 --- a/mfiles/tffilter.m +++ /dev/null @@ -1,97 +0,0 @@ -function y = tffilter(tfr,x,t,trace); -%TFFILTER Time frequency filtering of a signal. -% Y=TFFILTER(TFR,X,T,TRACE) filters the signal X -% with a non stationary filter. -% -% -% X : input signal (must be analytic). -% T : time instant(s) (default : 1:length(X)). -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% TFR : Wigner-Ville distribution of the filter -% frequency axis is graduated from 0.0 to 0.5. -% -% Example 1: -% Nt=128; t=1:Nt; sig=fmlin(Nt,0.05,0.3)+fmlin(Nt,0.2,0.45); -% sig(Nt/2)=sig(Nt/2)+8; figure(1);tfrwv(sig,t); -% Nf=128;freqs=0.5*(0:Nf-1).'/Nf; -% for tloop=1:Nt, -% rate=0.2*(tloop-1)/(Nt-1); -% H(:,tloop)=(0+rateNbPoints), - error('the values of T must be between 1 and xrow.'); -end; - -if trace, disp('non stationary signal filtering'); end; - -tfr=ifft(tfr); -y=zeros(tcol,1); - -for icol=1:tcol, - if trace, disprog(icol,tcol,10); end; - ti=t(icol); - valuestj=max([1 2-ti ti-N/2+1]):min([NbPoints 2*NbPoints-ti ti+N/2]); - %fprintf('%f %f %f\n',ti, min(valuestj),max(valuestj)); - - for tj=valuestj, - tmid = fix(0.5*(ti+tj)); - tdiff= fix(ti-tj); indice= rem(N+tdiff,N)+1; - %fprintf('%g %g %g\n',tj,tmid,indice); - y(icol,1)=y(icol,1)+tfr(indice,tmid)*x(tj); - end; -end; - -if trace, fprintf('\n'); end; diff --git a/mfiles/tfrbj.m b/mfiles/tfrbj.m deleted file mode 100644 index e3cb5b7..0000000 --- a/mfiles/tfrbj.m +++ /dev/null @@ -1,123 +0,0 @@ -function [tfr,t,f] = tfrbj(x,t,N,g,h,trace); -%TFRBJ Born-Jordan time-frequency distribution. -% [TFR,T,F]=TFRBJ(X,T,N,G,H,TRACE) computes the Born-Jordan -% distribution of a discrete-time signal X, or the -% cross Born-Jordan representation between two signals. -% -% X : signal if auto-BJ, or [X1,X2] if cross-BJ. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G : time smoothing window, G(0) being forced to 1. -% (default : Hamming(N/10)). -% H : frequency smoothing window, H(0) being forced to 1. -% (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, TFRBJ runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example : -% sig=fmlin(128,0.05,0.3)+fmlin(128,0.15,0.4); -% g=tftb_window(9,'Kaiser'); h=tftb_window(27,'Kaiser'); -% t=1:128; tfrbj(sig,t,128,g,h,1); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, May-August 1994, July 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -if (nargin <= 2), - N=xrow; -elseif (N<0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -hlength=floor(N/4); hlength=hlength+1-rem(hlength,2); -glength=floor(N/10); glength=glength+1-rem(glength,2); - -if (nargin == 1), - t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 2)|(nargin == 3), - g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 4), - h = tftb_window(hlength); trace = 0; -elseif (nargin == 5), - trace = 0; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('T must only have one row'); -end; - -[grow,gcol]=size(g); Lg=(grow-1)/2; % g=g/sum(g); -if (gcol~=1)|(rem(grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1); -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('Born-Jordan distribution'); end; -for icol=1:tcol, - ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]); - if trace, disprog(icol,tcol,10); end; - tfr(1,icol)= x(ti,1) .* conj(x(ti,xcol)); - - for tau=1:taumax, - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - g2=g(Lg+1+points); g2=g2/sum(g2); - R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol))); - tfr( 1+tau,icol)=h(Lh+tau+1)*R; - R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))); - tfr(N+1-tau,icol)=h(Lh-tau+1)*R; - end; - - tau=round(N/2); - if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh), - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - g2=g(Lg+1+points); g2=g2/sum(g2); - tfr(tau+1,icol) = 0.5 * ... - (h(Lh+tau+1)*sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)))+... - h(Lh-tau+1)*sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol)))); - end; -end; - -if trace, fprintf('\n'); end; -tfr= fft(tfr); -if (xcol==1), tfr=real(tfr); end ; - -if (nargout==0), - tfrqview(tfr,x,t,'tfrbj',g,h); -elseif (nargout==3), - f=(0.5*(0:N-1)/N)'; -end; diff --git a/mfiles/tfrbud.m b/mfiles/tfrbud.m deleted file mode 100644 index 4c6fea9..0000000 --- a/mfiles/tfrbud.m +++ /dev/null @@ -1,137 +0,0 @@ -function [tfr,t,f] = tfrbud(x,t,N,g,h,sigma,trace); -%TFRBUD Butterworth time-frequency distribution. -% [TFR,T,F]=TFRBUD(X,T,N,G,H,SIGMA,TRACE) computes the Butterworth -% distribution of a discrete-time signal X, or the -% cross Butterworth representation between two signals. -% -% X : signal if auto-BUD, or [X1,X2] if cross-BUD. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G : time smoothing window, G(0) being forced to 1. -% (default : Hamming(N/4)). -% H : frequency smoothing window, H(0) being forced to 1. -% (default : Hamming(N/4)). -% SIGMA : kernel width (default : 1). -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% TFR : time-frequency representation. When called without -% output arguments, TFRBUD runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example : -% sig=fmlin(128,0.05,0.3)+fmlin(128,0.15,0.4); -% g=tftb_window(9,'Kaiser'); h=tftb_window(27,'Kaiser'); -% t=1:128; tfrbud(sig,t,128,g,h,3.6,1); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, May-August 1994, July 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -if (nargin <= 2), - N=xrow; -elseif (N<0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -hlength=floor(N/4); hlength=hlength+1-rem(hlength,2); -glength=floor(N/10);glength=glength+1-rem(glength,2); - -if (nargin == 1), - t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); sigma = 1.0; trace = 0; -elseif (nargin == 2)|(nargin == 3), - g = tftb_window(glength); h = tftb_window(hlength); sigma = 1.0; trace = 0; -elseif (nargin == 4), - h = tftb_window(hlength); sigma = 1.0; trace = 0; -elseif (nargin == 5), - sigma = 1.0; trace = 0; -elseif (nargin == 6), - trace = 0; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('t must only have one row'); -end; - -[grow,gcol]=size(g); Lg=(grow-1)/2; -if (gcol~=1)|(rem(grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1); -if (hcol~=1)|(rem(hrow,2)==0), - error('H must be a smoothing window with odd length'); -end; - -if (sigma<=0.0), - error('SIGMA must be strictly positive'); -end; - -taumax = min([round(N/2),Lh]); tau = 1:taumax; points = -Lg:Lg; -BudKer = exp(-kron( abs(points.'), 1.0 ./ (2.0*tau/sqrt(sigma)))); -BudKer = diag(g) * BudKer; - -tfr= zeros (N,tcol) ; -if trace, disp('Butterworth distribution'); end; -for icol=1:tcol, - ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]); - if trace, disprog(icol,tcol,10); end; - tfr(1,icol)= x(ti,1) .* conj(x(ti,xcol)); - - for tau=1:taumax, - points= -min([Lg,xrow-ti-tau]):min([Lg,ti-tau-1]); - g2 = BudKer(Lg+1+points,tau); g2=g2/sum(g2); - R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol))); - tfr( 1+tau,icol)=h(Lh+tau+1)*R; - R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))); - tfr(N+1-tau,icol)=h(Lh-tau+1)*R; - end; - - tau=round(N/2); - if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh), - points= -min([Lg,xrow-ti-tau]):min([Lg,ti-tau-1]); - g2 = BudKer(Lg+1+points,tau); g2=g2/sum(g2); - tfr(tau+1,icol) = 0.5 * ... - (h(Lh+tau+1)*sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)))+... - h(Lh-tau+1)*sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol)))); - end; -end; - -clear BudKer; - -if trace, fprintf('\n'); end; -tfr= fft(tfr); -if (xcol==1), tfr=real(tfr); end ; - -if (nargout==0), - tfrqview(tfr,x,t,'tfrbud',g,h,sigma); -elseif (nargout==3), - f=(0.5*(0:N-1)/N)'; -end; - diff --git a/mfiles/tfrcw.m b/mfiles/tfrcw.m deleted file mode 100644 index 0511c70..0000000 --- a/mfiles/tfrcw.m +++ /dev/null @@ -1,141 +0,0 @@ -function [tfr,t,f] = tfrcw(x,t,N,g,h,sigma,trace); -%TFRCW Choi-Williams time-frequency distribution. -% [TFR,T,F]=TFRCW(X,T,N,G,H,SIGMA,TRACE) computes the Choi-Williams -% distribution of a discrete-time signal X, or the -% cross Choi-Williams representation between two signals. -% -% X : signal if auto-CW, or [X1,X2] if cross-CW. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G : time smoothing window, G(0) being forced to 1. -% (default : Hamming(N/10)). -% H : frequency smoothing window, H(0) being forced to 1. -% (default : Hamming(N/4)). -% SIGMA : kernel width (default : 1) -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% TFR : time-frequency representation. When called without -% output arguments, TFRCW runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example: -% sig=fmlin(128,0.05,0.3)+fmlin(128,0.15,0.4); -% g=tftb_window(9,'Kaiser'); h=tftb_window(27,'Kaiser'); -% t=1:128; tfrcw(sig,t,128,g,h,3.6,1); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, May-August 1994, July 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -if (nargin <= 2), - N=xrow; -elseif (N<0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -hlength=floor(N/4); hlength=hlength+1-rem(hlength,2); -glength=floor(N/10);glength=glength+1-rem(glength,2); - -if (nargin == 1), - t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); sigma = 1.0; trace = 0; -elseif (nargin == 2)|(nargin == 3), - g = tftb_window(glength); h = tftb_window(hlength); sigma = 1.0; trace = 0; -elseif (nargin == 4), - h = tftb_window(hlength); sigma = 1.0; trace = 0; -elseif (nargin == 5), - sigma = 1.0; trace = 0; -elseif (nargin == 6), - trace = 0; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('T must only have one row'); -end; - -[grow,gcol]=size(g); Lg=(grow-1)/2; -if (gcol~=1)|(rem(grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1); -if (hcol~=1)|(rem(hrow,2)==0), - error('H must be a smoothing window with odd length'); -end; - -if (sigma<=0.0), - error('SIGMA must be strictly positive'); -end; - -normfac = 16.0*pi/sigma; spreadfac = 16.0/sigma; -taumax = min([round(N/2),Lh]); tau = 1:taumax; points = -Lg:Lg; -CWKer = exp(-kron( points.' .^2, 1.0 ./ (spreadfac*tau.^2))); -CWKer = diag(g) * CWKer; - -tfr= zeros (N,tcol) ; -if trace, disp('Choi-Williams distribution'); end; -for icol=1:tcol, - ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]); - if trace, disprog(icol,tcol,10); end; - tfr(1,icol)= x(ti,1) .* conj(x(ti,xcol)); - - for tau=1:taumax, - points= -min([Lg,xrow-ti-tau]):min([Lg,ti-tau-1]); - g2 = CWKer(Lg+1+points,tau); g2=g2/sum(g2); - R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol))); - tfr( 1+tau,icol)=h(Lh+tau+1)*R; - R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))); - tfr(N+1-tau,icol)=h(Lh-tau+1)*R; - end; - - tau=round(N/2); - if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh), - points= -min([Lg,xrow-ti-tau]):min([Lg,ti-tau-1]); - g2 = CWKer(Lg+1+points,tau); g2=g2/sum(g2); - tfr(tau+1,icol) = 0.5 * ... - (h(Lh+tau+1)*sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)))+... - h(Lh-tau+1)*sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol)))); - end; -end; - -clear CWKer; - -if trace, fprintf('\n'); end; -tfr= fft(tfr); -if (xcol==1), tfr=real(tfr); end ; - -if (nargout==0), - tfrqview(tfr,x,t,'tfrcw',g,h,sigma); -elseif (nargout==3), - f=(0.5*(0:N-1)/N)'; -end; - - - - diff --git a/mfiles/tfrgrd.m b/mfiles/tfrgrd.m deleted file mode 100644 index 09074ec..0000000 --- a/mfiles/tfrgrd.m +++ /dev/null @@ -1,149 +0,0 @@ -function [tfr,t,f] = tfrgrd(x,t,N,g,h,rs,MoverN,trace); -%TFRGRD Generalized rectangular time-frequency distribution. -% [TFR,T,F]=TFRGRD(X,T,N,G,H,RS,MOVERN,TRACE) computes the -% Generalized Rectangular distribution of a discrete-time signal X, -% or the cross GRD representation between two signals. -% -% X : signal if auto-GRD, or [X1,X2] if cross-GRD. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G : time smoothing window, G(0) being forced to 1. -% (default : Hamming(N/10)). -% H : frequency smoothing window, H(0) being forced to 1. -% (default : Hamming(N/4)). -% RS : kernel width (default : 1). -% MOVERN : dissymmetry ratio (default : 1). -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% TFR : time-frequency representation. When called without -% output arguments, TFRGRD runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example : -% sig=fmlin(128,0.05,0.3)+fmlin(128,0.15,0.4); -% g=tftb_window(9,'Kaiser'); h=tftb_window(27,'Kaiser'); -% t=1:128; tfrgrd(sig,t,128,g,h,36,1/5,1); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, May-August 1994, July 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('at least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -if (nargin <= 2), - N=xrow; -elseif (N<0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -hlength=floor(N/4); hlength=hlength+1-rem(hlength,2); -glength=floor(N/10);glength=glength+1-rem(glength,2); - -if (nargin == 1), - t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); - rs = 1.0; MoverN = 1.0; trace = 0; -elseif (nargin == 2)|(nargin == 3), - g = tftb_window(glength); h = tftb_window(hlength); - rs = 1.0; MoverN = 1.0; trace = 0; -elseif (nargin == 4), - h = tftb_window(hlength); rs = 1.0; MoverN = 1.0; trace = 0; -elseif (nargin == 5), - rs = 1.0; MoverN = 1.0; trace = 0; -elseif (nargin == 6), - MoverN = 1.0; trace = 0; -elseif (nargin == 7), - trace = 0; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('T must only have one row'); -end; - -[grow,gcol]=size(g); Lg=(grow-1)/2; -if (gcol~=1)|(rem(grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1); -if (hcol~=1)|(rem(hrow,2)==0), - error('H must be a smoothing window with odd length'); -end; - -if (rs<=0.0), - error('RS must be strictly positive'); -end; - -if (MoverN<=0.0), - error('MOVERN must be strictly positive'); -end; - -taumax = min([round(N/2)-1,Lh]); tau = 1:taumax; points = -Lg:Lg; -GRDKer = sinc(-kron( 2.0*rs*points.', 1.0 ./ (2.0*tau).^MoverN )); -GRDKer = diag(g) * GRDKer; - -tfr= zeros (N,tcol) ; -if trace, disp('Generalized rectangular distribution'); end; -for icol=1:tcol, - ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]); - if trace, disprog(icol,tcol,10); end; - tfr(1,icol)= x(ti,1) .* conj(x(ti,xcol)); - - for tau=1:taumax, - points= -min([Lg,xrow-ti-tau]):min([Lg,ti-tau-1]); - g2 = GRDKer(Lg+1+points,tau) ; g2=g2/sum(g2); - R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol))); - tfr( 1+tau,icol)=h(Lh+tau+1)*R; - R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))); - tfr(N+1-tau,icol)=h(Lh-tau+1)*R; - end; - - tau=round(N/2); - if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh), - points= -min([Lg,xrow-ti-tau]):min([Lg,ti-tau-1]); - g2 = GRDKer(Lg+1+points,tau); g2=g2/sum(g2); - tfr(tau+1,icol) = 0.5 * ... - (h(Lh+tau+1)*sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)))+... - h(Lh-tau+1)*sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol)))); - end; -end; - -clear GRDKer; - -if trace, fprintf('\n'); end; -tfr= fft(tfr); -if (xcol==1), tfr=real(tfr); end ; - -if (nargout==0), - tfrqview(tfr,x,t,'tfrgrd',g,h,rs,MoverN); -elseif (nargout==3), - f=(0.5*(0:N-1)/N)'; -end; - - - - diff --git a/mfiles/tfristft.m b/mfiles/tfristft.m deleted file mode 100644 index 55a085d..0000000 --- a/mfiles/tfristft.m +++ /dev/null @@ -1,85 +0,0 @@ -function [x,t] = tfristft(tfr,t,h,trace); -%TFRISTFT Inverse Short time Fourier transform. -% [X,T]=TFRSTFT(tfr,T,H,TRACE) computes the inverse short-time -% Fourier transform of a discrete-time signal X. This function -% may be used for time-frequency synthesis of signals. -% -% X : signal. -% T : time instant(s) (default : 1:length(X)). -% H : frequency smoothing 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 decomposition (complex values). The -% frequency axis is graduated from -0.5 to 0.5. -% -% Example : -% t=200+(-128:127); sig=[fmconst(200,0.2);fmconst(200,0.4)]; -% h=hamming(57); tfr=tfrstft(sig,t,256,h,1); -% sigsyn=tfristft(tfr,t,h,1); -% plot(t,abs(sigsyn-sig(t))) -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, November 1996. -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin<3), - error('At least 3 parameters required'); -elseif (nargin==3), - trace=0; -end; - -[N,NbPoints]=size(tfr); -[trow,tcol] =size(t); -[hrow,hcol] =size(h); Lh=(hrow-1)/2; - -if (trow~=1), - error('T must only have one row'); -elseif (hcol~=1)|(rem(hrow,2)==0), - error('H must be a smoothing window with odd length'); -elseif (2^nextpow2(N)~=N), - fprintf('For a faster computation, N should be a power of two\n'); -elseif (tcol~=NbPoints) - error('tfr should have as many columns as t has rows.'); -end; - -Deltat=t(2:tcol)-t(1:tcol-1); -Mini=min(Deltat); Maxi=max(Deltat); -if (Mini~=1) & (Maxi~=1), - error('The tfr must be computed at each time sample.'); -end; - -h=h/norm(h); - -if trace, disp('Inverse Short-time Fourier transform'); end; -tfr=ifft(tfr); - -x=zeros(tcol,1); - -for icol=1:tcol, - if trace, disprog(icol,tcol,10); end; - valuestj=max([1,icol-N/2,icol-Lh]):min([tcol,icol+N/2,icol+Lh]); - for tj=valuestj, - tau=icol-tj; indices= rem(N+tau,N)+1; - % fprintf('%g %g %g\n',tj,tau,indices); - x(icol,1)=x(icol,1)+tfr(indices,tj)*h(Lh+1+tau); - end; - x(icol,1)=x(icol,1)/sum(abs(h(Lh+1+icol-valuestj)).^2); -end; - -if trace, fprintf('\n'); end; diff --git a/mfiles/tfrmhs.m b/mfiles/tfrmhs.m deleted file mode 100644 index 3b1552a..0000000 --- a/mfiles/tfrmhs.m +++ /dev/null @@ -1,115 +0,0 @@ -function [tfr,t,f] = tfrmhs(x,t,N,g,h,trace); -%TFRMHS Margenau-Hill-Spectrogram time-frequency distribution. -% [TFR,T,F]=TFRMHS(X,T,N,G,H,TRACE) computes the Margenau-Hill-Spectrogram -% distribution of a discrete-time signal X, or the cross -% Margenau-Hill-Spectrogram representation between two signals. -% -% X : Signal if auto-MHS, or [X1,X2] if cross-MHS. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G,H : analysis windows, normalized so that the representation -% preserves the signal energy. -% (default : Hamming(N/10) and Hamming(N/4)). -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% TFR : time-frequency representation. When called without -% output arguments, TFRMHS runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example: -% sig=fmlin(128,0.1,0.4); g=tftb_window(21,'Kaiser'); -% h=tftb_window(63,'Kaiser'); tfrmhs(sig,1:128,64,g,h,1); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, May-August 1994, July 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -if (nargin <= 2), - N=xrow; -elseif (N<0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -hlength=floor(N/4); hlength=hlength+1-rem(hlength,2); -glength=floor(N/10);glength=glength+1-rem(glength,2); - -if (nargin == 1), - t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 2)|(nargin == 3), - g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 4), - h = tftb_window(hlength); trace = 0; -elseif (nargin == 5), - trace = 0; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('T must only have one row'); -end; - -[grow,gcol]=size(g); Lg=(grow-1)/2; -if (gcol~=1)|(rem(grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1); -if (hcol~=1)|(rem(hrow,2)==0), - error('H must be a smoothing window with odd length'); -end; - -Lgh=min(Lg,Lh); points=-Lgh:Lgh; -Kgh=sum(h(Lh+1+points).*conj(g(Lg+1+points))); h=h/Kgh; - -tfr= zeros (N,tcol); tfr2= zeros(N,tcol); -if trace, disp('Pseudo Margenau-Hill distribution'); end; -for icol=1:tcol, - ti= t(icol); - if trace, disprog(icol,tcol,10); end; - tau=-min([round(N/2)-1,Lg,ti-1]):min([round(N/2)-1,Lg,xrow-ti]); - indices= rem(N+tau,N)+1; - tfr(indices,icol)=x(ti+tau,1).*conj(g(Lg+1+tau)); - tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]); - indices= rem(N+tau,N)+1; - tfr2(indices,icol)=x(ti+tau,xcol).*conj(h(Lh+1+tau)); -end; -if trace, fprintf('\n'); end; -tfr=real(fft(tfr).*conj(fft(tfr2))); - -if (nargout==0), - tfrqview(tfr,x,t,'tfrmhs',g,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; - - diff --git a/mfiles/tfrmmce.m b/mfiles/tfrmmce.m deleted file mode 100644 index de2aae8..0000000 --- a/mfiles/tfrmmce.m +++ /dev/null @@ -1,96 +0,0 @@ -function [tfr,t,f] = tfrmmce(x,h,t,N,trace); -%TFRMMCE Minimum mean cross-entropy combination of spectrograms. -% [TFR,T,F]=TFRMMCE(X,H,T,N,TRACE) computes the minimum mean -% cross-entropy combination of spectrograms using as -% windows the columns of the matrix H. -% -% X : signal. -% H : frequency smoothing windows, the H(:,i) being normalized -% so as to be of unit energy. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)) -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% TFR : time-frequency representation. When called without -% output arguments, TFRMMCE runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example : -% sig=fmlin(128,0.1,0.4); h=zeros(19,3); -% h(10+(-5:5),1)=tftb_window(11); h(10+(-7:7),2)=tftb_window(15); -% h(10+(-9:9),3)=tftb_window(19); tfrmmce(sig,h); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, August 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -[xrow,xcol] = size(x); -if (nargin < 2), - error('At least 2 parameters are required'); -elseif (nargin==2), - t=1:xrow; N=xrow; trace=0; -elseif (nargin==3), - N=xrow; trace=0; -elseif (nargin==4), - trace=0; -end; -[trow,tcol] = size(t); - -if (xcol~=1), - error('X must have only one column'); -elseif (trow~=1), - error('T must only have one row'); -elseif (N<=0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N & nargin==5), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; -if (rem(hrow,2)==0), - error('H must have an odd number of lines'); -elseif hcol==1, - error('H must have at least two columns'); -end; -h=h*diag(1.0 ./ sqrt(sum(abs(h).^2))); - -tfr= zeros (N,tcol); -slides= zeros(N,hcol); -if trace, disp('MMCE 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; slides= zeros(N,hcol); - if trace, disprog(icol,tcol,10); end; - for ih=1:hcol, - slides(indices,ih)=x(ti+tau).*conj(h(Lh+1+tau,ih)); - end; - slides=abs(fft(slides)).^2; - tfr(:,icol)=prod(slides')' .^(1/hcol); -end; - -tfr=tfr*(sum(abs(x(t)).^2)/sum(sum(tfr))); -if trace, fprintf('\n'); end; - -if (nargout==0), - tfrqview(tfr,x,t,'tfrmmce',h); -elseif (nargout==3), - f=fftshift((-(N/2):(N/2)-1)/N)'; -end; - diff --git a/mfiles/tfrparam.m b/mfiles/tfrparam.m deleted file mode 100644 index a49a0c5..0000000 --- a/mfiles/tfrparam.m +++ /dev/null @@ -1,97 +0,0 @@ -function [tfr,t,f] = tfrparam(x,t,N,p,L,trace,Rxtype,method,q); -% [tfr,t,f]=tfrparam(x,t,N,p,L,trace,Rxtype,method,q) parametric time-frequency -% representation. -% -% X : signal -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : max(256,length(X)) ). -% p : last autocorrelation lag (default : 2 ). -% L : length of the window around the analyzed time sample. Must be odd. -% (default : max(51,round(xrow/10)) ). -% Rxtype : choice of the correlation matrix algorithm (default : 'fbhermitian') -% possible values are 'hermitian', 'fbhermitian', 'burg' or 'fbburg' -% method : can be either 'ar', 'periodogram', 'capon', 'capnorm', 'lagunas', -% or 'genlag'. -% q : parameter for the generalized Lagunas method. -% -% example: -% -% x=fmconst(256,0.1)+fmlin(256,0.15,0.30)+fmlin(256,0.2,0.45); -% figure(1); tfrparam(x,1:256,256,3,21,1,'fbhermitian','ar'); -% figure(2); tfrparam(x,1:256,256,15,61,1,'fbhermitian','periodogram'); -% figure(3); tfrparam(x,1:256,256,3,21,1,'fbhermitian','capon'); -% figure(4); tfrparam(x,1:256,256,3,21,1,'fbhermitian','lagunas'); - -% F. Auger, july, november 1998, april 99. -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (nargin==1), - t=1:xrow; N=max(256,xrow); p=2; L=max(51,round(xrow/10)); trace=0; Rxtype='fbhermitian'; method='ar'; -elseif (nargin==2), - N=max(256,xrow); p=2; L=max(51,round(xrow/10)); trace=0; Rxtype='fbhermitian'; method='ar'; -elseif (nargin==3), - p=2; L=max(51,round(xrow/10)); trace=0; Rxtype='fbhermitian'; method='ar'; -elseif (nargin==4), - L=max(51,round(xrow/10)); trace=0; Rxtype='fbhermitian'; method='ar'; -elseif (nargin==5), - trace=0; Rxtype='fbhermitian'; method='ar'; -elseif (nargin==6), - Rxtype='fbhermitian'; method='ar'; -elseif (nargin==7), - method='ar'; -end; - -[trow,tcol] = size(t); - -if xcol>1, - error('x must be a column vector'); -elseif p<1, - error('p must be greater than 0'); -end; - -if (rem(L,2)==0), % L must be odd - L=L+1; fprintf('L corrected to L+1\n'); -end; -Lhalf=(L-1)/2; - -tfr= zeros (N,tcol) ; -if trace, disp('Tfrparam'); end; -for icol=1:tcol, - ti= t(icol); - if ti-Lhalf >= 1, timin=ti-Lhalf; else timin=1; end; - if ti+Lhalf <= xrow, timax=ti+Lhalf; else timax=xrow; end; - if trace, disprog(icol,tcol,10); end; - Realp=min (p,timax-timin); - Rx=correlmx(x(timin:timax),Realp,Rxtype); - if strcmp(upper(method),'GENLAG'), - [tfr(:,icol),f]=parafrep(Rx,N,method,q); - else - [tfr(:,icol),f]=parafrep(Rx,N,method); - end -end; -if trace, fprintf('\n'); end; - -if (nargout==0), - imagesc(t,f,10.0*log10(tfr)); axis('xy'); - title(['tfrparam, ' method, ... - ', p=', num2str(p), ... - ', L=', num2str(L), ', ', ... - Rxtype, ' autocorrelation matrix']); drawnow; -end; - - -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA diff --git a/mfiles/tfrrgab.m b/mfiles/tfrrgab.m deleted file mode 100644 index b8dd1cb..0000000 --- a/mfiles/tfrrgab.m +++ /dev/null @@ -1,172 +0,0 @@ -function [tfr,rtfr,hat] = tfrrgab(x,t,N,Nh,trace,K); -%TFRRGAB Reassigned Gabor spectrogram time-frequency distribution. -% [TFR,RTFR,HAT] = TFRRGAB(X,T,N,NH,TRACE,K) -% computes the Gabor spectrogram and its reassigned version. -% This particular window (a Gaussian window) allows a 20 % faster -% algorithm than the TFRRSP function. -% -% X : analysed signal -% T : the time instant(s) (default : 1:length(X)) -% N : number of frequency bins (default : length(X)) -% NH : length of the gaussian window (default : N/4)) -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% K : value at both extremities (default 0.001) -% TFR, : time-frequency representation and its reassigned -% RTFR version. When called without output arguments, -% TFRRGAB runs TFRQVIEW. -% HAT : Complex matrix of the reassignment vectors. -% -% Example : -% sig=fmlin(128,0.1,0.4); tfrrgab(sig,1:128,128,19,1); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, May-July 1994, July 1995. -% Copyright (c) 1996 by CNRS(France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (nargin <= 2), - N=xrow; -end; - -hlength=floor(N/4); -hlength=hlength+1-rem(hlength,2); - -if (nargin == 1), - t=1:xrow; -end; - -if (nargin <= 3), - Nh=hlength; trace=0; K=0.001; -elseif (nargin == 4), - trace = 0; K=0.001; -elseif (nargin == 5), - K= 0.001; -end; - -if (N<0), - error('N must be greater than zero'); -end; -[trow,tcol] = size(t); -if (xcol~=1), - error('X must have only one column'); -elseif (trow~=1), - error('T must only have one row'); -elseif (2^nextpow2(N)~=N & nargin==6), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -if (rem(Nh,2)==0), - error('Nh must be odd'); -elseif length(Nh)~=1, - error('Nh must be a scalar'); -end; - -Nh2=Nh-2; -TFTBcontinue=1; -while TFTBcontinue, - Nh2=Nh2+2; - h=tftb_window(Nh2,'gauss',K^((Nh2-1)^2 /(Nh-1)^2)); - TFTBcontinue=(h(Nh2)*(Nh2-1)>2*K); -end; - -K=K^((Nh2-1)^2 /(Nh-1)^2); Nh=Nh2; Lh=(Nh-1)/2; -h=h; Th=h.*[-Lh:Lh]'; - -if (tcol==1), - Dt=1; -else - Deltat=t(2:tcol)-t(1:tcol-1); - Mini=min(Deltat); Maxi=max(Deltat); - if (Mini~=Maxi), - error('The time instants must be regularly sampled.'); - else - Dt=Mini; - end; - clear Deltat Mini Maxi; -end; - -tfr= zeros(N,tcol); tf2= zeros(N,tcol); tf3= zeros(N,tcol); -if trace, disp('Gabor spectrogram'); end; - -for icol=1:tcol, - if trace, disprog(icol,tcol,10); end; - 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; - norm_h=norm(h(Lh+1+tau)); - tfr(indices,icol)=x(ti+tau).*conj( h(Lh+1+tau)) /norm_h; - tf2(indices,icol)=x(ti+tau).*conj(Th(Lh+1+tau)) /norm_h; -end ; -tfr=fft(tfr); tf2=fft(tf2); -tfr=tfr(:); tf2=tf2(:); tf3=tf3(:); -avoid_warn=find(tfr~=0.0); -tf3(avoid_warn)=round(imag(2*log(K)*N*tf2(avoid_warn)./tfr(avoid_warn)/(2.0*pi*Lh^2))); -tf2(avoid_warn)=round(real(tf2(avoid_warn)./tfr(avoid_warn)/Dt)); -tfr=abs(tfr).^2; -if trace, fprintf ('\nreassignment: \n'); end; -tfr=reshape(tfr,N,tcol); -tf2=reshape(tf2,N,tcol); -tf3=reshape(tf3,N,tcol); - -rtfr= zeros(N,tcol); -Ex=mean(abs(x(min(t):max(t))).^2); Threshold=1.0e-6*Ex; -for icol=1:tcol, - if trace, disprog(icol,tcol,10); end; - for jcol=1:N, - if abs(tfr(jcol,icol))>Threshold, - icolhat= icol + tf2(jcol,icol); - icolhat=min(max(icolhat,1),tcol); - jcolhat= jcol - tf3(jcol,icol); - %while (jcolhat<1),jcolhat=jcolhat+N; end; - %while (jcolhat>N),jcolhat=jcolhat-N; end; - jcolhat=rem(rem(jcolhat-1,N)+N,N)+1; - rtfr(jcolhat,icolhat)=rtfr(jcolhat,icolhat) + tfr(jcol,icol) ; - tf2(jcol,icol)=jcolhat + j * icolhat; - else - tf2(jcol,icol)=inf*(1+j); - rtfr(jcol,icol)=rtfr(jcol,icol) + tfr(jcol,icol) ; - end; - end; -end; - -if trace, fprintf('\n'); end; -clear tf3; -if (nargout==0), - TFTBcontinue=1; - while (TFTBcontinue==1), - choice=menu ('Choose the representation:',... - 'stop',... - 'Gabor spectrogram',... - 'reassigned Gabor spectrogram'); - if (choice==1), TFTBcontinue=0; - elseif (choice==2), - Q=round(tcol*N/xrow); - tfrqview(tfr,x,t,'tfrgabor',tcol,Q,h); - elseif (choice==3), - tfrqview(rtfr,x,t,'tfrrgab',Nh); - end; - end; -elseif (nargout>2), - hat=tf2; -end; - diff --git a/mfiles/tfrri.m b/mfiles/tfrri.m deleted file mode 100644 index e4152ce..0000000 --- a/mfiles/tfrri.m +++ /dev/null @@ -1,87 +0,0 @@ -function [tfr,t,f] = tfrri(x,t,N,trace); -%TFRRI Rihaczek time-frequency distribution. -% [TFR,T,F]=TFRRI(X,T,N,TRACE) computes the Rihaczek distribution -% of a discrete-time signal X, -% or the cross Rihaczek representation between two signals. -% -% X : signal if auto-Ri, or [X1,X2] if cross-Ri. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% TFR : time-frequency representation. When called without -% output arguments, TFRRI applies TFRQVIEW on the real -% part of the distribution, which is equal to the -% Margenau-Hill distribution. -% F : vector of normalized frequencies. -% -% Example : -% sig=fmlin(128,0.1,0.4); tfrri(sig); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, May-August 1994, July 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; - -[xrow,xcol] = size(x); -if (nargin == 1), - t=1:xrow; N=xrow; trace=0; -elseif (nargin == 2), - N=xrow; trace=0; -elseif (nargin == 3), - 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 & nargin==4), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -tfr= zeros (N,tcol) ; -if trace, disp('Rihaczek distribution'); end; -for icol=1:tcol, - ti= t(icol); tau=-min([N-ti,xrow-ti]):(ti-1); - indices= rem(N+tau,N)+1; - if trace, disprog(icol,tcol,10); end; - tfr(indices,icol)=x(ti,1)*conj(x(ti-tau,xcol)); -end; -if trace, fprintf('\n'); end; -tfr= fft(tfr); - -if (nargout==0), - fprintf('if your goal is a graphical output, you should rather use tfrmh.\n'); - tfrqview(real(tfr),x,t,'tfrmh'); -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; - diff --git a/mfiles/tfrridb.m b/mfiles/tfrridb.m deleted file mode 100644 index 881d938..0000000 --- a/mfiles/tfrridb.m +++ /dev/null @@ -1,130 +0,0 @@ -function [tfr,t,f] = tfrridb(x,t,N,g,h,trace); -%TFRRIDB Reduced Interference Distribution with Bessel kernel. -% [TFR,T,F]=TFRRIDB(X,T,N,G,H,TRACE) Reduced Interference -% Distribution with a kernel based on the Bessel function of the -% first kind. -% TFRRIDB computes either the distribution of a discrete-time -% signal X, or the cross representation between two signals. -% -% X : signal if auto-RIDB, or [X1,X2] if cross-RIDB. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G : time smoothing window, G(0) being forced to 1. -% (default : Hamming(N/10)). -% H : frequency smoothing window, H(0) being forced to 1. -% (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, TFRRIDB runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example : -% sig=[fmlin(128,0.05,0.3)+fmlin(128,0.15,0.4);zeros(64,1)]; -% g=tftb_window(31,'rect'); h=tftb_window(63,'rect'); -% for t=120:1:130; t,plot(tfrridb(sig,t,128,g,h,0)); end; -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, August 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -if (nargin <= 2), - N=xrow; -elseif (N<0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N & nargin==6), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -hlength=floor(N/4); hlength=hlength+1-rem(hlength,2); -glength=floor(N/10);glength=glength+1-rem(glength,2); - -if (nargin == 1), - t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 2)|(nargin == 3), - g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 4), - h = tftb_window(hlength); trace = 0; -elseif (nargin == 5), - trace = 0; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('T must only have one row'); -end; - -[grow,gcol]=size(g); Lg=(grow-1)/2; g=g/g(Lg+1); -if (gcol~=1)|(rem(grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1); -if (hcol~=1)|(rem(hrow,2)==0), - error('H must be a smoothing window with odd length'); -end; - -taumax = min([round(N/2)-1,Lh]); - -tfr= zeros (N,tcol) ; -if trace, disp('Guo-Durand-Lee distribution'); end; -for icol=1:tcol, - ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]); - if trace, disprog(icol,tcol,10); end; - tfr(1,icol)= g(Lg+1) * x(ti,1) .* conj(x(ti,xcol)); - - for tau=1:taumax, - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - Ker= sqrt(1.0 - (points/tau).^2)'; - g2 = g(Lg+1+points) .* Ker; - if (sum(g2)>eps), g2=g2/sum(g2); end; - R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol))); - tfr( 1+tau,icol)=h(Lh+tau+1)*R; - R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))); - tfr(N+1-tau,icol)=h(Lh-tau+1)*R; - end; - - tau=round(N/2); - if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh), - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - Ker= sqrt(1.0 - (points/tau).^2)'; - g2 = g(Lg+1+points) .* Ker; g2=g2/sum(g2); - tfr(tau+1,icol) = 0.5 * ... - (h(Lh+tau+1)*sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)))+... - h(Lh-tau+1)*sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol)))); - end; -end; - -if trace, fprintf('\n'); end; -tfr= fft(tfr); -if (xcol==1), tfr=real(tfr); end ; - -if (nargout==0), - tfrqview(tfr,x,t,'tfrridb',g,h); -elseif (nargout==3), - f=(0.5*(0:N-1)/N)'; -end; diff --git a/mfiles/tfrridbn.m b/mfiles/tfrridbn.m deleted file mode 100644 index 6f14e52..0000000 --- a/mfiles/tfrridbn.m +++ /dev/null @@ -1,132 +0,0 @@ -function [tfr,t,f] = tfrridbn(x,t,N,g,h,trace); -%TFRRIDBN Reduced Interference Distribution with a binomial kernel. -% [TFR,T,F]=TFRRIDBN(X,T,N,G,H,TRACE) Reduced Interference -% Distribution with a kernel based on the binomial coefficients. -% TFRRIDBN computes either the distribution of a discrete-time -% signal X, or the cross representation between two signals. -% -% X : signal if auto-RIDBN, or [X1,X2] if cross-RIDBN. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G : time smoothing window, G(0) being forced to 1. -% (default : Hamming(N/10)). -% H : frequency smoothing window, H(0) being forced to 1. -% (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, TFRRIDBN runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example : -% sig=[fmlin(128,.05,.3)+fmlin(128,.15,.4)]; tfrridbn(sig); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, June 1996. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -if (nargin <= 2), - N=xrow; -elseif (N<0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N & nargin==6), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -hlength=floor(N/4); hlength=hlength+1-rem(hlength,2); -glength=floor(N/10);glength=glength+1-rem(glength,2); - -if (nargin == 1), - t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 2)|(nargin == 3), - g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 4), - h = tftb_window(hlength); trace = 0; -elseif (nargin == 5), - trace = 0; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('T must only have one row'); -end; - -[grow,gcol]=size(g); Lg=(grow-1)/2; g=g/g(Lg+1); -if (gcol~=1)|(rem(grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1); -if (hcol~=1)|(rem(hrow,2)==0), - error('H must be a smoothing window with odd length'); -end; - -taumax = min([round(N/2)-1,Lh]); - -tfr= zeros (N,tcol) ; -if trace, disp('Binomial RID distribution'); end; -for icol=1:tcol, - ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]); - if trace, disprog(icol,tcol,10); end; - tfr(1,icol)= g(Lg+1) * x(ti,1) .* conj(x(ti,xcol)); - - Ker=[1]; - for tau=1:taumax, - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - Ker= ([Ker; 0; 0]+2*[0; Ker; 0]+[0; 0; Ker])/4.0; - g2 = g(Lg+1+points) .* Ker(tau+points+1); - if (sum(g2)>eps), g2=g2/sum(g2); end; - R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol))); - tfr( 1+tau,icol)=h(Lh+tau+1)*R; - R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))); - tfr(N+1-tau,icol)=h(Lh-tau+1)*R; - end; - - tau=round(N/2); - if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh), - Ker=ones(2*tau+1,1); - for p=1:2*tau, - Ker(p+1)=Ker(p)*(2*tau-p+1)/p; - end; - Ker=Ker/sum(Ker); - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - g2 = g(Lg+1+points) .* Ker(tau+points+1); g2=g2/sum(g2); - tfr(tau+1,icol) = 0.5 * ... - (h(Lh+tau+1)*sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)))+... - h(Lh-tau+1)*sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol)))); - end; -end; - -if trace, fprintf('\n'); end; -tfr= fft(tfr); -if (xcol==1), tfr=real(tfr); end ; - -if (nargout==0), - tfrqview(tfr,x,t,'tfrridbn',g,h); -elseif (nargout==3), - f=(0.5*(0:N-1)/N)'; -end; diff --git a/mfiles/tfrridh.m b/mfiles/tfrridh.m deleted file mode 100644 index 7e10651..0000000 --- a/mfiles/tfrridh.m +++ /dev/null @@ -1,130 +0,0 @@ -function [tfr,t,f] = tfrridh(x,t,N,g,h,trace); -%TFRRIDH Reduced Interference Distribution with Hanning kernel. -% [TFR,T,F]=TFRRIDH(X,T,N,G,H,TRACE) Reduced Interference -% Distribution with a kernel based on the Hanning window. -% TFRRIDH computes either the distribution of a discrete-time -% signal X, or the cross representation between two signals. -% -% X : signal if auto-RIDH, or [X1,X2] if cross-RIDH. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G : time smoothing window, G(0) being forced to 1. -% (default : Hamming(N/10)). -% H : frequency smoothing window, H(0) being forced to 1. -% (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, TFRRIDH runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example : -% sig=[fmlin(128,0.05,0.3)+fmlin(128,0.15,0.4);zeros(64,1)]; -% g=tftb_window(31,'rect'); h=tftb_window(63,'rect'); -% for t=120:1:130; t,plot(tfrridh(sig,t,128,g,h,0)); end; -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, May-August 1994, August 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -if (nargin <= 2), - N=xrow; -elseif (N<0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N & nargin==6), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -hlength=floor(N/4); hlength=hlength+1-rem(hlength,2); -glength=floor(N/10);glength=glength+1-rem(glength,2); - -if (nargin == 1), - t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 2)|(nargin == 3), - g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 4), - h = tftb_window(hlength); trace = 0; -elseif (nargin == 5), - trace = 0; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('T must only have one row'); -end; - -[grow,gcol]=size(g); Lg=(grow-1)/2; g=g/g(Lg+1); -if (gcol~=1)|(rem(grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1); -if (hcol~=1)|(rem(hrow,2)==0), - error('H must be a smoothing window with odd length'); -end; - -taumax = min([round(N/2)-1,Lh]); - -tfr= zeros (N,tcol) ; -if trace, disp('Hamming RID distribution'); end; -for icol=1:tcol, - ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]); - if trace, disprog(icol,tcol,10); end; - tfr(1,icol)= g(Lg+1) * x(ti,1) .* conj(x(ti,xcol)); - - for tau=1:taumax, - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - Ker= (1.0+cos(pi*points/tau))'; - g2 = g(Lg+1+points) .* Ker; - if (sum(g2)>eps), g2=g2/sum(g2); end; - R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol))); - tfr( 1+tau,icol)=h(Lh+tau+1)*R; - R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))); - tfr(N+1-tau,icol)=h(Lh-tau+1)*R; - end; - - tau=round(N/2); - if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh), - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - Ker= (1.0+cos(pi*points/tau))'; - g2 = g(Lg+1+points) .* Ker; g2=g2/sum(g2); - tfr(tau+1,icol) = 0.5 * ... - (h(Lh+tau+1)*sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)))+... - h(Lh-tau+1)*sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol)))); - end; -end; - -if trace, fprintf('\n'); end; -tfr= fft(tfr); -if (xcol==1), tfr=real(tfr); end ; - -if (nargout==0), - tfrqview(tfr,x,t,'tfrridh',g,h); -elseif (nargout==3), - f=(0.5*(0:N-1)/N)'; -end; - diff --git a/mfiles/tfrridt.m b/mfiles/tfrridt.m deleted file mode 100644 index f1d0b6d..0000000 --- a/mfiles/tfrridt.m +++ /dev/null @@ -1,131 +0,0 @@ -function [tfr,t,f] = tfrridt(x,t,N,g,h,trace); -%TFRRIDT Reduced Interference Distribution with triangular kernel. -% [TFR,T,F]=TFRRIDT(X,T,N,G,H,TRACE) Reduced Interference -% Distribution with a kernel based on the triangular (or Bartlett) -% window. -% TFRRIDT computes either the distribution of a discrete-time -% signal X, or the cross representation between two signals. -% -% X : signal if auto-RIDT, or [X1,X2] if cross-RIDT. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G : time smoothing window, G(0) being forced to 1. -% (default : Hamming(N/10)). -% H : frequency smoothing window, H(0) being forced to 1. -% (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, TFRRIDT runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example : -% sig=[fmlin(128,0.05,0.3)+fmlin(128,0.15,0.4);zeros(64,1)]; -% g=tftb_window(31,'rect'); h=tftb_window(63,'rect'); -% for t=120:1:130; t,plot(tfrridt(sig,t,128,g,h,0)); end; -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, May-August 1994, August 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -if (nargin <= 2), - N=xrow; -elseif (N<0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N & nargin==6), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -hlength=floor(N/4); hlength=hlength+1-rem(hlength,2); -glength=floor(N/10);glength=glength+1-rem(glength,2); - -if (nargin == 1), - t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 2)|(nargin == 3), - g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 4), - h = tftb_window(hlength); trace = 0; -elseif (nargin == 5), - trace = 0; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('T must only have one row'); -end; - -[grow,gcol]=size(g); Lg=(grow-1)/2; g=g/g(Lg+1); -if (gcol~=1)|(rem(grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1); -if (hcol~=1)|(rem(hrow,2)==0), - error('H must be a smoothing window with odd length'); -end; - -taumax = min([round(N/2)-1,Lh]); - -tfr= zeros (N,tcol) ; -if trace, disp('triangular RID distribution'); end; -for icol=1:tcol, - ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]); - if trace, disprog(icol,tcol,10); end; - tfr(1,icol)= g(Lg+1) * x(ti,1) .* conj(x(ti,xcol)); - - for tau=1:taumax, - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - Ker= (1.0-abs(points/tau))'; - g2 = g(Lg+1+points) .* Ker; - if (sum(g2)>eps), g2=g2/sum(g2); end; - R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol))); - tfr( 1+tau,icol)=h(Lh+tau+1)*R; - R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))); - tfr(N+1-tau,icol)=h(Lh-tau+1)*R; - end; - - tau=round(N/2); - if (ti<=xrow-tau)&(ti>=tau+1)&(tau<=Lh), - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - Ker= (1.0-abs(points/tau))'; - g2 = g(Lg+1+points) .* Ker; g2=g2/sum(g2); - tfr(tau+1,icol) = 0.5 * ... - (h(Lh+tau+1)*sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol)))+... - h(Lh-tau+1)*sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol)))); - end; -end; - -if trace, fprintf('\n'); end; -tfr= fft(tfr); -if (xcol==1), tfr=real(tfr); end ; - -if (nargout==0), - tfrqview(tfr,x,t,'tfrridt',g,h); -elseif (nargout==3), - f=(0.5*(0:N-1)/N)'; -end; - diff --git a/mfiles/tfrrstan.m b/mfiles/tfrrstan.m deleted file mode 100644 index a49bdef..0000000 --- a/mfiles/tfrrstan.m +++ /dev/null @@ -1,164 +0,0 @@ -function [stan,rtfr,hat] = tfrrstan(x,t,N,G,h,trace); -%TFRRSTAN Reassigned Stankovic distribution. -% [TFR,RTFR,HAT] = TFRRSTAN(X,T,N,H,TRACE) -% computes the Stankovic distribution and its reassigned version. -% -% X : analysed signal. -% T : the time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G : frequency averaging window -% h : stft window (default : Hamming(N/4)). -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% TFR, : time-frequency representation and its reassigned -% RTFR version. When called without output arguments, -% TFRRSTAN runs TFRQVIEW. -% HAT : Complex matrix of the reassignment vectors. -% -% Example : -% sig=fmlin(128,0.1,0.4); t=1:2:128; -% G=tftb_window(9,'hanning'); h=tftb_window(61,'hanning'); tfrrstan(sig,t,128,G,h,1); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, August, September 1997. -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (nargin <= 2), - N=xrow; -end; - -hlength=floor(N/4); -hlength=hlength+1-rem(hlength,2); - -if (nargin == 1), - t=1:xrow; G=[0.25; 0.5; 0.25]; h = tftb_window(hlength); trace=0; -elseif (nargin == 2)|(nargin == 3), - G=[0.25; 0.5; 0.25]; h = tftb_window(hlength); trace=0; -elseif (nargin == 4), - h = tftb_window(hlength); trace = 0; -elseif (nargin == 5), - trace = 0; -end; - -[Grow,Gcol]=size(G); LG=(Grow-1)/2; -if (Gcol~=1)|(rem(Grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -if (N<0), - error('N must be greater than zero'); -end; -[trow,tcol] = size(t); -if (xcol~=1), - error('X must have only one column'); -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; - -if (tcol==1), - Dt=1; -else - Deltat=t(2:tcol)-t(1:tcol-1); - Mini=min(Deltat); Maxi=max(Deltat); - if (Mini~=Maxi), - error('The time instants must be regularly sampled.'); - else - Dt=Mini; - end; - clear Deltat Mini Maxi; -end; - -stan = zeros(N,tcol); -rtfr = zeros(N,tcol); -hat = zeros(N,tcol); - -if trace, disp('Stankovic distribution (with reassignement)'); end; -Ex=mean(abs(x(min(t):max(t))).^2); Threshold=1.0e-3*Ex; -Dh=dwindow(h); Th=h.*[-Lh:Lh]'; -for icol=1:tcol, - if trace, disprog(icol,tcol,10); end; - 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; - tfr= zeros(N,3); - tfr(indices,1)=x(ti+tau).*conj( h(Lh+1-tau)); - tfr(indices,2)=x(ti+tau).*conj(Th(Lh+1-tau)); - tfr(indices,3)=x(ti+tau).*conj(Dh(Lh+1-tau)); - tfr=fft(tfr); - stan(:,icol)=G(LG+1) * abs(tfr(:,1)).^2; - - for jcol=1:N, - stanTh=G(LG+1)*tfr(jcol,1)*conj(tfr(jcol,2)); - stanDh=G(LG+1)*tfr(jcol,1)*conj(tfr(jcol,3)); - for kstan=1:min(N/2-1,LG), - stanbefore=rem(rem(jcol-kstan-1,N)+N,N)+1; - stanafter =rem(rem(jcol+kstan-1,N)+N,N)+1; - stan(jcol,icol)= stan(jcol,icol) ... - + G(LG+1-kstan)*tfr(stanbefore,1)*conj(tfr(stanafter ,1)) ... - + G(LG+1+kstan)*tfr(stanafter ,1)*conj(tfr(stanbefore,1)); - stanTh= stanTh + G(LG+1-kstan)*tfr(stanbefore,1)*conj(tfr(stanafter ,2)) ... - + G(LG+1+kstan)*tfr(stanafter ,1)*conj(tfr(stanbefore,2)); - stanDh= stanDh + G(LG+1-kstan)*tfr(stanbefore,1)*conj(tfr(stanafter ,3)) ... - + G(LG+1+kstan)*tfr(stanafter ,1)*conj(tfr(stanbefore,3)); - end; - stan(jcol,icol)=real(stan(jcol,icol)); - if abs(stan(jcol,icol))>Threshold, - icolhat = round(icol - real(stanTh/stan(jcol,icol))/Dt); - jcolhat = round(jcol - N*imag(stanDh/stan(jcol,icol)/(2.0*pi))); - - jcolhat= rem(rem(jcolhat-1,N)+N,N)+1; - icolhat= min(max(icolhat,1),tcol); - - rtfr(jcolhat,icolhat)=rtfr(jcolhat,icolhat) + stan(jcol,icol) ; - hat(jcol,icol)= jcolhat + j * icolhat; - %fprintf('%12.3f %12.3f , %12.3f %12.3f \n',jcol,icol,jcolhat,icolhat); - else - rtfr(jcol,icol)=rtfr(jcol,icol) + stan(jcol,icol); - hat(jcol,icol)= jcol + j * icol; - end; - end; - -end ; - -if trace, fprintf('\n'); end; - -if (nargout==0), - TFTBcontinue=1; - while (TFTBcontinue==1), - choice=menu ('Choose the representation:',... - 'stop',... - 'Stankovic distribution',... - 'reassigned Stankovic distribution'); - if (choice==1), TFTBcontinue=0; - elseif (choice==2), - tfrqview(stan,x,t,'type1'); - elseif (choice==3), - tfrqview(rtfr,x,t,'type1'); - end; - end; -end; diff --git a/mfiles/tfrsave.m b/mfiles/tfrsave.m deleted file mode 100644 index 9b6aa56..0000000 --- a/mfiles/tfrsave.m +++ /dev/null @@ -1,103 +0,0 @@ -function tfrsave(name,tfr,method,sig,t,f,p1,p2,p3,p4,p5); -%TFRSAVE Save the parameters of a time-frequency representation. -% TFRSAVE(NAME,TFR,METHOD,SIG,T,F,P1,P2,P3,P4,P5) saves the -% parameters of a time-frequency representation in the file -% NAME.mat. Two additional parameters are saved : TfrQView and -% TfrView. If you load the file 'name.mat' and do eval(TfrQView), you -% will restart the display session under tfrqview ; if you do -% eval(TfrView), you will display the representation by means of -% tfrview. -% -% NAME : name of the mat-file (less than 8 characters). -% TFR : time-frequency representation (MxN). -% METHOD : chosen representation. -% SIG : signal from which the TFR was obtained -% T : time instant(s) (default : (1:N)). -% F : frequency bins (default : linspace(0,0.5,M)). -% P1..P5 : optional parameters : run the file TFRPARAM(METHOD) -% to know the meaning of P1..P5 for your method. -% -% Example : -% sig=fmlin(64); tfr=tfrwv(sig); -% tfrsave('wigner',tfr,'TFRWV',sig,1:64); -% clear; load wigner; eval(TfrQView); -% -% See also TFRQVIEW, TFRVIEW, TFRPARAM. - -% F. Auger, August 1994 - O. Lemoine, June 1996. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -[M,N]=size(tfr); - -if (nargin < 3), - error('At least 3 parameters are required'); -elseif nargin==3, - sig=[]; t=1:N; f=(0.5*(0:M-1)/M); -elseif nargin==4, - t=1:N; f=(0.5*(0:M-1)/M); -elseif nargin==5, - f=(0.5*(0:M-1)/M); -end; - -method=upper(method); -namedflt=name; - -while (length(name)>8), - disp('The name must have less than 8 characters'); - namedflt=name(1:8); - nameStr=[' Name of the MAT file [',namedflt,'] : ']; - name=input(nameStr,'s'); -end - -if name==[], - name=namedflt; -end - -if (nargin==3), - t=1:N; - TfrQView=['tfrqview(tfr,[],t,method)']; - TfrView =['clf;tfrview(tfr,sig,t,method,param,map)']; - eval(['save ',name,' tfr t f method TfrQView TfrView']); -elseif (nargin==4), - TfrQView=['tfrqview(tfr,[],t,method)']; - TfrView =['clf;tfrview(tfr,sig,t,method,param,map)']; - eval(['save ',name,' tfr t f method TfrQView TfrView']); -elseif (nargin==5), - TfrQView=['tfrqview(tfr,sig,t,method)']; - TfrView =['clf;tfrview(tfr,sig,t,method,param,map)']; - eval(['save ',name,' tfr sig t f method TfrQView TfrView']); -elseif (nargin==6), - TfrQView=['tfrqview(tfr,sig,t,method,p1)']; - TfrView =['clf; tfrview(tfr,sig,t,method,param,map,p1)']; - eval(['save ',name,' tfr sig t f method p1 TfrQView TfrView']); -elseif (nargin==7), - TfrQView=['tfrqview(tfr,sig,t,method,p1,p2)']; - TfrView =['clf; tfrview(tfr,sig,t,method,param,map,p1,p2)']; - eval(['save ',name,' tfr sig t f method p1 p2 TfrQView TfrView']); -elseif (nargin==8), - TfrQView=['tfrqview(tfr,sig,t,method,p1,p2,p3)']; - TfrView =['clf; tfrview(tfr,sig,t,method,param,map,p1,p2,p3)']; - eval(['save ',name,' tfr sig t f method p1 p2 p3 TfrQView TfrView']); -elseif (nargin==9), - TfrQView=['tfrqview(tfr,sig,t,method,p1,p2,p3,p4)']; - TfrView =['clf; tfrview(tfr,sig,t,method,param,map,p1,p2,p3,p4)']; - eval(['save ',name,' tfr sig t f method p1 p2 p3 p4 TfrQView TfrView']); -elseif (nargin==10), - TfrQView=['tfrqview(tfr,sig,t,method,p1,p2,p3,p4,p5)']; - TfrView =['clf; tfrview(tfr,sig,t,method,param,map,p1,p2,p3,p4,p5)']; - eval(['save ',name,' tfr sig t f method p1 p2 p3 p4 p5 TfrQView TfrView']); -end; diff --git a/mfiles/tfrscalt.m b/mfiles/tfrscalt.m deleted file mode 100644 index 6add82f..0000000 --- a/mfiles/tfrscalt.m +++ /dev/null @@ -1,135 +0,0 @@ -%function tfrscalt -%TFRSCALT Unit test for the time-frequency representation TFRSCALO. - -% O. Lemoine - June 1996. - -% We test each property of the corresponding TFR : - -N=128; - -% Covariance by translation in time -t1=60; t2=70; f=0.3; W=0; -sig1=amgauss(N,t1).*fmconst(N,f,t1); -sig2=amgauss(N,t2).*fmconst(N,f,t2); -tfr1=tfrscalo(sig1,1:N,W,0.1,0.4,128); -tfr2=tfrscalo(sig2,1:N,W,0.1,0.4,128); -[tr,tc]=size(tfr1); -nu=round(f*(tc-1)*2)+1; -tfr=tfr1-tfr2(:,modulo((1:tc)-t1+t2,tc)); -if any(any(abs(tfr)>sqrt(eps))), - error('tfrscalo test 1 failed'); -end - - -% Covariance by dilation -t=N/2; f=0.2; T=2*sqrt(N); a=2; W=8; -sig1=amgauss(N,t,T).*fmconst(N,f,t); -sig2=amgauss(a*N,a*t,T*a).*fmconst(a*N,f/a,a*t); -[tfr1,t1,f1]=tfrscalo(sig1,1:N ,W,0.01,0.49,N); -[tfr2,t2,f2]=tfrscalo(sig2,1:a*N,W,0.01,0.49,N); -Max1=max(max(tfr1)); Max2=max(max(tfr2)); -[I1,J1]=find(tfr1==Max1); [I2,J2]=find(tfr2==Max2); -if abs(f1(I1)-a*f2(I2))>1e-2 | J2~=a*J1, - error('tfrscalo test 2 failed'); -end - - -% Reality of the TFR -sig=noisecg(N); W=5; -tfr=tfrscalo(sig,1:N,W,0.01,0.5,N); -if sum(any(abs(imag(tfr))>sqrt(eps)))~=0, - error('tfrscalo test 3 failed'); -end - - -% Energy conservation -sig=fmsin(N,.1,.4); W=6; Nf=2*N ; -[tfr,t,f]=tfrscalo(sig,1:N,W,0.01,0.49,Nf); -Es=norm(sig)^2/Nf; -Etfr=integ2d(tfr,t,f)/N; -if abs(Es-Etfr)>sqrt(eps), - error('tfrscalo test 4 failed'); -end - - -% Positivity -if any(any(tfr<0)), - error('tfrscalo test 5 failed'); -end - - -% Same energy in the time-scale plane for 2 gaussian atoms at different scales -sig=amgauss(256).*(fmconst(256,.15)+fmconst(256,.35)); -[tfr,t,f]=tfrscalo(sig,1:256,12,.01,.49,512); -int1=integ2d(tfr(1:430,:),t,f(1:430)); -int2=integ2d(tfr(431:512,:),t,f(431:512)); -if abs(int1-int2)>1e-4, - error('tfrscalo test 6 failed'); -end - - -N=127; - -% Covariance by dilation -t=ceil(N/2); f=0.2; T=2*sqrt(N); a=2; W=8; -sig1=amgauss(N,t,T).*fmconst(N,f,t); -sig2=amgauss(a*N,a*t,T*a).*fmconst(a*N,f/a,a*t); -[tfr1,t1,f1]=tfrscalo(sig1,1:N ,W,0.01,0.49,N); -[tfr2,t2,f2]=tfrscalo(sig2,1:a*N,W,0.01,0.49,N); -Max1=max(max(tfr1)); Max2=max(max(tfr2)); -[I1,J1]=find(tfr1==Max1); [I2,J2]=find(tfr2==Max2); -if abs(f1(I1(1))-a*f2(I2))>1e-2 | J2~=a*J1(1), - error('tfrscalo test 7 failed'); -end - - -% Reality of the TFR -sig=noisecg(N); W=5; -tfr=tfrscalo(sig,1:N,W,0.01,0.5,N); -if sum(any(abs(imag(tfr))>sqrt(eps)))~=0, - error('tfrscalo test 8 failed'); -end - - -% Energy conservation -sig=fmsin(N,.1,.4); W=6; Nf=2*N+1 ; -[tfr,t,f]=tfrscalo(sig,1:N,W,0.01,0.49,Nf); -SP = fft(hilbert(real(sig))); -indmin = 1+round(0.01*(N-2)); -indmax = 1+round(0.49*(N-2)); -SPana = SP(indmin:indmax); -Es=SPana'*SPana/Nf; -Etfr=integ2d(tfr,t,f); -if abs(Es-Etfr)>1e-2, - error('tfrscalo test 9 failed'); -end - - -% Positivity -if any(any(tfr<0)), - error('tfrscalo test 10 failed'); -end - - -% Same energy in the time-scale plane for 2 gaussian atoms at different scales -sig=amgauss(N).*(fmconst(N,.15)+fmconst(N,.35)); -[tfr,t,f]=tfrscalo(sig,1:N,12,.01,.49,2*N+1); -int1=integ2d(tfr(1:210,:),t,f(1:210)); -int2=integ2d(tfr(211:255,:),t,f(211:255)); -if abs(int1-int2)>5e-4, - error('tfrscalo test 11 failed'); -end - -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA diff --git a/mfiles/tfrspbk.m b/mfiles/tfrspbk.m deleted file mode 100644 index 5d2cc30..0000000 --- a/mfiles/tfrspbk.m +++ /dev/null @@ -1,272 +0,0 @@ -function [tfr,t,f]=tfrspbk(X,time,K,nh0,ng0,fmin,fmax,N,trace); -%TFRSPBK Smoothed Pseudo K-Bertrand time-frequency distribution. -% [TFR,T,F]=TFRSPBK(X,T,K,NH0,NG0,FMIN,FMAX,N,TRACE) -% generates the auto- or cross- Smoothed Pseudo K-Bertrand -% distribution. -% -% X : signal (in time) to be analyzed. If X=[X1 X2], TFRSPBK -% computes the cross-Smoothed Pseudo K-Bertrand distribution. -% (Nx=length(X)). -% T : time instant(s) on which the TFR is evaluated. TIME must -% be a uniformly sampled vector whose elements are between 1 -% and Nx. (default : 1:Nx). -% K : label of the K-Bertrand distribution. The distribution with -% parametrization function -% lambdak(u,K) = (K (exp(-u)-1)/(exp(-Ku)-1))^(1/(K-1)) -% is computed (default : 0). -% K=-1 : Smoothed pseudo (active) Unterberger distribution -% K=0 : Smoothed pseudo Bertrand distribution -% K=1/2: Smoothed pseudo D-Flandrin distribution -% K=2 : Affine smoothed pseudo Wigner-Ville distribution. -% NH0 : half length of the analyzing wavelet at coarsest scale. -% A Morlet wavelet is used. NH0 controles the frequency -% smoothing of the smoothed pseudo K-Bertrand distribution. -% (default : sqrt(Nx)). -% NG0 : half length of the time smoothing window. -% NG0 = 0 corresponds to the Pseudo K-Bertrand distribution. -% (default : 0). -% FMIN,FMAX : respectively lower and upper frequency bounds of -% the analyzed signal. These parameters fix the equivalent -% frequency bandwidth (expressed in Hz). When unspecified, you -% have to enter them at the command line from the plot of the -% spectrum. FMIN and FMAX must be >0 and <=0.5. -% N : number of analyzed voices (default : Nx). -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% TFR : time-frequency matrix containing the coefficients of the -% decomposition (abscissa correspond to uniformly sampled time, -% and ordinates correspond to a geometrically sampled -% frequency). First row of TFR corresponds to the lowest -% frequency. When called without output arguments, TFRSPBK -% runs TFRQVIEW. -% F : vector of normalized frequencies (geometrically sampled -% from FMIN to FMAX). -% -% Example : -% sig=altes(64,0.1,0.45); tfrspbk(sig); -% -% See also TFRBERT, TFRUNTAC, TFRUNTPA, TFRSCALO, TFRDFLA, TFRASPW. - -% P. Goncalves, October 95 - O. Lemoine, June 1996. -% Copyright (c) 1995 Rice University -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least one parameter required'); -end; - -[xrow,xcol] = size(X); -if nargin<=8, trace=0; end - -if (nargin == 1), - time=1:xrow; K=0; nh0=sqrt(xrow); ng0=0; -elseif (nargin == 2), - K=0; nh0=sqrt(xrow); ng0=0; -elseif (nargin == 3), - nh0=sqrt(xrow); ng0=0; -elseif (nargin == 4), - ng0=0; -elseif (nargin == 6), - disp('FMIN will not be taken into account. Determine it with FMAX'); - disp(' from the following plot of the spectrum.'); -elseif (nargin == 7), - N=xrow; -end; - -[trow,tcol] = size(time); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -elseif (trow~=1), - error('TIME must only have one row'); -end; - -Mt=length(X); - -if trace, disp('Smoothed Pseudo K-Bertrand distribution'); end; - -if xcol==1, - X1=X; - X2=X; -else - X1=X(:,1); - X2=X(:,2); -end -s1 = real(X1); -s2 = real(X2); - -if rem(Mt,2)~=0, - s1 = [s1;0]; - s2 = [s2;0]; - M = (Mt+1)/2; -else - M = Mt/2; -end ; - -t = [-nh0:nh0-1]; -Tmin = 1 ; -Tmax = 2*nh0 ; -T = Tmax-Tmin ; - -if nargin<=6, % fmin,fmax,N unspecified - STF1 = fft(fftshift(s1(min(time):max(time)))); Nstf=length(STF1); - sp1 = (abs(STF1(1:Nstf/2))).^2; Maxsp1=max(sp1); - STF2 = fft(fftshift(s2(min(time):max(time)))); - sp2 = (abs(STF2(1:Nstf/2))).^2; Maxsp2=max(sp2); - f = linspace(0,0.5,Nstf/2+1) ; f=f(1:Nstf/2); - plot(f,sp1) ; grid; hold on ; plot(f,sp2) ; hold off - xlabel('Normalized frequency'); - title('Analyzed signal energy spectrum'); - axis([0 1/2 0 1.2*max(Maxsp1,Maxsp2)]) ; - indmin=min(find(sp1>Maxsp1/100)); - indmax=max(find(sp1>Maxsp1/100)); - fmindflt=max([0.01 0.05*fix(f(indmin)/0.05)]); - fmaxdflt=0.05*ceil(f(indmax)/0.05); - txtmin=['Lower frequency bound [',num2str(fmindflt),'] : ']; - txtmax=['Upper frequency bound [',num2str(fmaxdflt),'] : ']; - fmin = input(txtmin); fmax = input(txtmax); - if fmin==[], fmin=fmindflt; end - if fmax==[], fmax=fmaxdflt; end -end - -if fmin >= fmax - error('FMAX must be greater or equal to FMIN'); -elseif fmin<=0.0 | fmin>0.5, - error('FMIN must be > 0 and <= 0.5'); -elseif fmax<=0.0 | fmax>0.5, - error('FMAX must be > 0 and <= 0.5'); -end - -B = fmax-fmin ; -R = B/((fmin+fmax)/2) ; -Qte = fmax/fmin ; -umax = log(Qte); -Teq = nh0/(fmax*umax); -if Teq<2*nh0, - M0 = round((2*nh0^2)/Teq-nh0)+1; - MU = nh0+M0; - T2 = 2*MU-1; -else - M0 = 0; - MU = nh0; - T2 = 2*MU-1; -end; - -if nargin<=6, - Nq= ceil((B*T2*(1+2/R)*log((1+R/2)/(1-R/2)))/2); - Nmin = Nq-rem(Nq,2); - Ndflt = 2^nextpow2(Nmin); - Ntxt=['Number of frequency samples (>=',num2str(Nmin),') [',num2str(Ndflt),'] : ']; - N = input(Ntxt); - if N==[], N=Ndflt; end -end - -fmin_s = num2str(fmin) ; fmax_s = num2str(fmax) ; N_s = num2str(N) ; -if trace, - disp(['Frequency runs from ',fmin_s,' to ',fmax_s,' with ',N_s,' points']); -end - - -k = 1:N; -q = (fmax/fmin)^(1/(N)); -a = (exp((k-1).*log(q))); % a is an increasing scale vector. -geo_f(k) = fmin*a ; % geo_f is a geometrical increasing - % frequency vector. - -% Morlet wavelet decomposition computation -t0 = 1; -t1 = Mt; -Mtr = Mt; -z1 = hilbert(s1).'; -matxte1 = zeros(N,Mt); -z2 = hilbert(s2).'; -matxte2 = zeros(N,Mt); -nu0 = geo_f(N); -for ptr=1:N, - nha = round(nh0*a(ptr)); - nua = nu0/a(ptr); - ha = exp(-(2*log(10)/(nh0*a(ptr))^2)*(-nha:nha).^2).*exp(-i*2*pi*nua*(-nha:nha)); - detail1 = conv(z1(t0:t1),fliplr(ha)); - matxte1(N-ptr+1,:) = detail1(nha+1:length(detail1)-nha); - detail2 = conv(z2(t0:t1),fliplr(ha)); - matxte2(N-ptr+1,:) = detail2(nha+1:length(detail2)-nha); - % first row of matxte corresponds to the lowest frequency. -end; - - -% Pseudo-Bertrand distribution computation -tfr=zeros(N,tcol); -umin = -umax; -u=linspace(umin,umax,2*MU+1); -U(MU+1) = 0; -k = 1:2*N; -beta(k) = -1/(2*log(q))+(k-1)./(2*N*log(q)); -for m = 1:2*MU+1, - l1(m,:) = exp(-(2*i*pi*beta+1/2).*log(lambdak(u(m),K))); -end -if ng0==0 - decay = 0 ; -elseif ng0~=0 - gamma0 = ng0*fmax ; - alpha = - log(0.01)/gamma0^2 ; - u0 = sqrt(-alpha*log(-0.01*sqrt(alpha/pi))/pi^2) ; - decay = -log(0.01)*(umax/u0)^2/log(10) ; -end -if decay==Inf - G = zeros(1,2*MU) ; - G(MU+1) = 1 ; -elseif decay==0, - G=ones(1,2*MU); -else - Nb=2*MU+1; - G = amgauss(Nb,(Nb+1)/2,(Nb-1)*sqrt(pi/(decay*log(10)))/2).'; - G = G(1:2*MU) ; -end -xx = exp(-[0:N-1].*log(q)); -xx = xx(ones(1,2*MU),:).*G(ones(1,N),:)'; -indi=1; -for ti = time, - if trace, disprog(ti-time(1)+1,time(tcol)-time(1)+1,10); end - S1 = zeros(1,2*N); - S1(1:N) = matxte1(:,ti).'; - Mellin1 = fftshift(ifft(S1.*exp([0:2*N-1].*log(q)))) ; - S2 = zeros(1,2*N); - S2(1:N) = matxte2(:,ti).'; - Mellin2 = fftshift(ifft(S2.*exp([0:2*N-1].*log(q)))) ; - waf = zeros(2*MU,N) ; - MX1 = l1.*Mellin1(ones(1,2*MU+1),:) ; - X1 = fft(MX1.'); - X1 = X1(1:N,:).' ; - MX2 = l1.*Mellin2(ones(1,2*MU+1),:) ; - X2 = fft(MX2.'); - X2 = X2(1:N,:).'; - waf = real(X1(1:2*MU,:).*conj(X2(2*MU+1:-1:2,:)).*xx) ; - tfr(:,indi) = sum(waf).'; % first row of tfr corresponds to - indi = indi+1; % the lowest frequency. -end; - -t = time; -f = geo_f'; -tfr = tfr./integ2d(tfr,t,f)*sum(s1.*conj(s2)) ; - -disp(' '); - -if (nargout==0), - tfrqview(tfr,X,t,'tfrspbk',K,nh0,ng0,N,f); -end; - - - - diff --git a/mfiles/tfrsurf.m b/mfiles/tfrsurf.m deleted file mode 100644 index d0515f9..0000000 --- a/mfiles/tfrsurf.m +++ /dev/null @@ -1,68 +0,0 @@ -function [tfr2,OrderedSurfaces]=tfrsurf(tfr,threshold,keep,trace); -% [tfr2,OrderedSurfaces]=tfrsurf(tfr,threshold,keep,trace); -% extract from a time-frequency representation the biggest energy dots -% TFR : time-frequency representation. -% THRESHOLD : the energy threshold, in % -% KEEP : number of dots to keep -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% -% example : -% -% N=256; -% sig=fmlin(N,0.1,0.3)+fmlin(N,0.3,0.4)+2*fmlin(N,0.05,0.2).*amgauss(N,190,70); -% tfr=tfrwv(sig,1:N,128); -% [tfr2,OrderedSurfaces]=tfrsurf(tfr,5,3,1); -% figure(1);tfrview(tfr,sig,1:N,'tfrwv',[2 1 5 10 128 2 1 5]) -% title('original tfr'); -% figure(2);tfrview(tfr2,sig,1:N,'tfrwv',[2 1 5 10 128 2 1 5]); -% title('modified tfr'); -% figure(3);semilogy(1:10,OrderedSurfaces(1:10),'-',1:10,OrderedSurfaces(1:10),'o'); -% title('number of points of the 10 biggest dots'); -% -% See also imextract. - -% F. Auger, oct 1999 -% Copyright (c) CNRS - France 1999. -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -if nargin==1, - threshold=5; keep=10; trace=0; -elseif nargin==2, - keep=10; trace=0; -elseif nargin==3, - trace=0; -end; - -[Nbrow,Nbcol]=size(tfr); -TheMax=max(max(tfr)); -[EnergyDots,NbDots]=imextrac((tfr>=threshold*TheMax*0.01),trace); -Surfaces=zeros(1,NbDots+1); -for i=0:NbDots, - Surfaces(i+1)=length(find(EnergyDots==i)); -end; -[OrderedSurfaces,Indices]=sort(Surfaces(2:NbDots+1)); -OrderedSurfaces=fliplr(OrderedSurfaces); -Indices=fliplr(Indices); - -Binary=zeros(Nbrow,Nbcol); -for i=1:keep, - DotIndice=find(EnergyDots==Indices(i)); - Binary(DotIndice)=1; -end; - -tfr2=tfr.*Binary; - - diff --git a/mfiles/tfrzam.m b/mfiles/tfrzam.m deleted file mode 100644 index e85b116..0000000 --- a/mfiles/tfrzam.m +++ /dev/null @@ -1,115 +0,0 @@ -function [tfr,t,f] = tfrzam(x,t,N,g,h,trace); -%TFRZAM Zao-Atlas-Marks time-frequency distribution. -% [TFR,T,F]=TFRZAM(X,T,N,G,H,TRACE) computes the Zao-Atlas-Marks -% distribution of a discrete-time signal X, or the -% cross Zao-Atlas-Marks representation between two signals. -% -% X : signal if auto-ZAM, or [X1,X2] if cross-ZAM. -% T : time instant(s) (default : 1:length(X)). -% N : number of frequency bins (default : length(X)). -% G : time smoothing window, G(0) being forced to 1. -% (default : Hamming(N/10)). -% H : frequency smoothing window, H(0) being forced to 1. -% (default : Hamming(N/4)). -% TRACE : if nonzero, the progression of the algorithm is shown -% (default : 0). -% TFR : time-frequency representation. When called without -% outpout arguments, TFRZAM runs TFRQVIEW. -% F : vector of normalized frequencies. -% -% Example : -% sig=fmlin(128,0.05,0.3)+fmlin(128,0.15,0.4); -% g=tftb_window(9,'Kaiser'); h=tftb_window(27,'Kaiser'); -% tfrzam(sig,1:128,128,g,h,1); -% -% See also all the time-frequency representations listed in -% the file CONTENTS (TFR*) - -% F. Auger, May-August 1994, July 1995. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -if (nargin == 0), - error('At least 1 parameter required'); -end; -[xrow,xcol] = size(x); -if (xcol==0)|(xcol>2), - error('X must have one or two columns'); -end - -if (nargin <= 2), - N=xrow; -elseif (N<0), - error('N must be greater than zero'); -elseif (2^nextpow2(N)~=N), - fprintf('For a faster computation, N should be a power of two\n'); -end; - -hlength=floor(N/4); hlength=hlength+1-rem(hlength,2); -glength=floor(N/10);glength=glength+1-rem(glength,2); - -if (nargin == 1), - t=1:xrow; g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 2)|(nargin == 3), - g = tftb_window(glength); h = tftb_window(hlength); trace = 0; -elseif (nargin == 4), - h = tftb_window(hlength); trace = 0; -elseif (nargin == 5), - trace = 0; -end; - -[trow,tcol] = size(t); -if (trow~=1), - error('T must only have one row'); -end; - -[grow,gcol]=size(g); Lg=(grow-1)/2; g=g/g(Lg+1); -if (gcol~=1)|(rem(grow,2)==0), - error('G must be a smoothing window with odd length'); -end; - -[hrow,hcol]=size(h); Lh=(hrow-1)/2; h=h/h(Lh+1); -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('Zao-Atlas-Marks Representation'); end; -for icol=1:tcol, - ti= t(icol); taumax=min([ti+Lg-1,xrow-ti+Lg,round(N/2)-1,Lh]); - if trace, disprog(icol,tcol,10); end; - tfr(1,icol)= g(Lg+1) * x(ti,1) .* conj(x(ti,xcol)); - - for tau=1:taumax, - points= -min([tau,Lg,xrow-ti-tau]):min([tau,Lg,ti-tau-1]); - g2=g(Lg+1+points); - R=sum(g2 .* x(ti+tau-points,1) .* conj(x(ti-tau-points,xcol))); - tfr( 1+tau,icol)=h(Lh+tau+1)*R; - R=sum(g2 .* x(ti-tau-points,1) .* conj(x(ti+tau-points,xcol))); - tfr(N+1-tau,icol)=h(Lh-tau+1)*R; - end; -end; - -if trace, fprintf('\n'); end; -tfr= fft(tfr); -if (xcol==1), tfr=real(tfr); end ; - -if (nargout==0), - tfrqview(tfr,x,t,'tfrzam',g,h); -elseif (nargout==3), - f=(0.5*(0:N-1)/N)'; -end; - diff --git a/mfiles/zak.m b/mfiles/zak.m deleted file mode 100644 index 184ebfe..0000000 --- a/mfiles/zak.m +++ /dev/null @@ -1,58 +0,0 @@ -function dzt=zak(sig,N,M) -%ZAK Zak transform. -% DZT=ZAK(SIG,N,M) computes the Zak transform of signal SIG. -% -% SIG : Signal to be analysed (length(X)=N1). -% N : number of Zak coefficients in time (N1 must be a multiple -% of N) (default : closest integer towards sqrt(N1)). -% M : number of Zak coefficients in frequency (N1 must be a -% multiple of M) (default : M=N1/N). -% DZT : Output matrix (N,M) containing the discrete Zak transform. -% -% Example : -% sig=fmlin(256); DZT=zak(sig); -% imagesc(DZT); -% -% See also IZAK, TFRGABOR. - -% O. Lemoine - February 1996. -% Copyright (c) 1996 by CNRS (France). -% -% This program is free software; you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation; either version 2 of the License, or -% (at your option) any later version. -% -% This program is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with this program; if not, write to the Free Software -% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -N1=length(sig); - -if N1<=2, - error('SIG must have at least 3 elements'); -elseif nargin<1, - error('The number of parameters must be at least one'); -elseif nargin==1, - [N,M]=divider(N1); -elseif nargin==2, - N=N1/M; -elseif M<=1|N<=1|M>=N1|N>=N1, - error('M and N must be between 2 and N1-1'); -elseif rem(N1,M)~=0|rem(N1,N)~=0, - error('M and N must be such that N1 is a multiple of M and N'); -end - -Msig=zeros(M,N); -dzt=zeros(N,M); - -Msig=reshape(sig,M,N); - -dzt=fft(Msig.')/sqrt(N); - - \ No newline at end of file diff --git a/tftb/__init__.py b/tftb/__init__.py index f15dc3f..9b41e1f 100644 --- a/tftb/__init__.py +++ b/tftb/__init__.py @@ -1,6 +1,6 @@ import sys -__version__ = '0.1.3' +__version__ = '0.1.4' try: __TFTB__SETUP__ diff --git a/tftb/processing/base.py b/tftb/processing/base.py index be7c6ba..79b6347 100644 --- a/tftb/processing/base.py +++ b/tftb/processing/base.py @@ -11,6 +11,7 @@ """ import numpy as np +from scipy.signal import hamming import matplotlib.pyplot as plt @@ -18,7 +19,7 @@ class BaseTFRepresentation(object): isaffine = False - def __init__(self, signal, **kwargs): + def __init__(self, signal, timestamps=None, n_fbins=None, fwindow=None, **kwargs): """Create a base time-frequency representation object. :param signal: Signal to be analyzed. @@ -27,28 +28,25 @@ def __init__(self, signal, **kwargs): :return: BaseTFRepresentation object :rtype: """ - if (signal.ndim == 2) and (1 in signal.shape): + if signal.ndim > 1: signal = signal.ravel() self.signal = signal - timestamps = kwargs.get('timestamps') if timestamps is None: timestamps = np.arange(signal.shape[0]) self.ts = self.timestamps = timestamps - n_fbins = kwargs.get('n_fbins') if n_fbins is None: n_fbins = signal.shape[0] self.n_fbins = n_fbins - fwindow = kwargs.get('fwindow') if fwindow is None: fwindow = self._make_window() self.fwindow = fwindow if self.n_fbins % 2 == 0: - freqs = np.hstack((np.arange(self.n_fbins / 2), - np.arange(-self.n_fbins / 2, 0))) + freqs = np.r_[np.arange(self.n_fbins / 2), + np.arange(-self.n_fbins / 2, 0)] else: - freqs = np.hstack((np.arange((self.n_fbins - 1) / 2), - np.arange(-(self.n_fbins - 1) / 2, 0))) - self.freqs = freqs.astype(float) / self.n_fbins + freqs = np.r_[np.arange((self.n_fbins - 1) / 2), + np.arange(-(self.n_fbins - 1) / 2, 0)] + self.freqs = freqs / self.n_fbins self.tfr = np.zeros((self.n_fbins, self.ts.shape[0]), dtype=complex) def _get_spectrum(self): @@ -68,12 +66,9 @@ def _make_window(self): :rtype: array-like """ - h = np.floor(self.n_fbins / 4.0) - h += 1 - np.remainder(h, 2) - fwindow = np.hamming(int(h)) - # No need to normalize the window - # fwindow = fwindow / np.linalg.norm(fwindow) - return fwindow + h = self.n_fbins // 4 + h += 1 - (h % 2) + return hamming(h) def _plot_tfr(self, ax, kind, extent, contour_x=None, contour_y=None, levels=None, show_tf=True, cmap=plt.cm.gray): @@ -196,17 +191,16 @@ def plot(self, ax=None, kind='cmap', show=True, default_annotation=True, self._annotate_signal(axTime) self._annotate_spectrum(axSpec) else: - if (ax is None) and (kind != "surf"): - fig = plt.figure() - ax = fig.add_subplot(111) + fig = plt.figure() if kind == "cmap": + ax = fig.add_subplot(111) ax.imshow(self.tfr, aspect='auto', cmap=cmap, origin='lower', extent=extent, **kwargs) elif kind == "surf": from mpl_toolkits.mplot3d import Axes3D fig = plt.figure() - ax = fig.gca(projection="3d") + ax = fig.add_subplot(projection="3d") x = np.arange(self.signal.shape[0]) y = np.linspace(0, 0.5, self.signal.shape[0]) X, Y = np.meshgrid(x, y) @@ -215,7 +209,7 @@ def plot(self, ax=None, kind='cmap', show=True, default_annotation=True, ax.set_zlabel("Amplitude") elif kind == "wireframe": from mpl_toolkits.mplot3d import Axes3D # NOQA - ax = fig.gca(projection="3d") + ax = fig.add_subplot(projection="3d") x = np.arange(self.signal.shape[0]) y = np.linspace(0, 0.5, self.signal.shape[0]) X, Y = np.meshgrid(x, y) @@ -223,6 +217,7 @@ def plot(self, ax=None, kind='cmap', show=True, default_annotation=True, rstride=3, cstride=3) else: t, f = np.meshgrid(self.ts, np.linspace(0, 0.5, self.tfr.shape[0])) + ax = fig.add_subplot(111) ax.contour(t, f, self.tfr, **kwargs) if default_annotation: grid = kwargs.get('grid', True) diff --git a/tftb/processing/cohen.py b/tftb/processing/cohen.py index a3778f2..a03204b 100644 --- a/tftb/processing/cohen.py +++ b/tftb/processing/cohen.py @@ -240,16 +240,16 @@ def smoothed_pseudo_wigner_ville(signal, timestamps=None, freq_bins=None, freq_bins = signal.shape[0] if fwindow is None: - winlength = np.floor(freq_bins / 4.0) - winlength = winlength + 1 - np.remainder(winlength, 2) + winlength = freq_bins // 4 + winlength = winlength + 1 - (winlength % 2) from scipy.signal import hamming fwindow = hamming(int(winlength)) elif fwindow.shape[0] % 2 == 0: raise ValueError('The smoothing fwindow must have an odd length.') if twindow is None: - timelength = np.floor(freq_bins / 10.0) - timelength += 1 - np.remainder(timelength, 2) + timelength = freq_bins // 10 + timelength += 1 - (timelength % 2) from scipy.signal import hamming twindow = hamming(int(timelength)) elif twindow.shape[0] % 2 == 0: @@ -261,7 +261,7 @@ def smoothed_pseudo_wigner_ville(signal, timestamps=None, freq_bins=None, for icol in range(timestamps.shape[0]): ti = timestamps[icol] taumax = min([ti + lg - 1, signal.shape[0] - ti + lg, - np.round(freq_bins / 2.0) - 1, lh]) + round(freq_bins / 2.0) - 1, lh]) points = np.arange(-min([lg, signal.shape[0] - ti]), min([lg, ti - 1]) + 1).astype(int) lg_points = (lg + points).astype(int) @@ -278,9 +278,9 @@ def smoothed_pseudo_wigner_ville(signal, timestamps=None, freq_bins=None, idx1 = (ti + tau - points - 1).astype(int) idx2 = (ti - tau - points - 1).astype(int) R = np.sum(g2 * signal[idx1] * np.conj(signal[idx2])) - tfr[1 + tau, icol] = fwindow[(lh + tau + 1).astype(int)] * R + tfr[1 + tau, icol] = fwindow[lh + tau + 1] * R R = np.sum(g2 * signal[idx2] * np.conj(signal[idx1])) - tfr[freq_bins - tau - 1, icol] = fwindow[(lh - tau + 1).astype(int)] * R + tfr[freq_bins - tau - 1, icol] = fwindow[lh - tau + 1] * R tau = np.round(freq_bins / 2.0) if (ti <= signal.shape[0] - tau) and (ti >= tau + 1) and (tau <= lh): points = np.arange(-min([lg, signal.shape[0] - ti - tau]), diff --git a/tftb/processing/linear.py b/tftb/processing/linear.py index dda3592..8615525 100644 --- a/tftb/processing/linear.py +++ b/tftb/processing/linear.py @@ -118,9 +118,9 @@ def gabor(signal, n_coeff=None, q_oversample=None, window=None): :rtype: tuple """ if n_coeff is None: - n_coeff = divider(signal.shape[0]) + n_coeff, _ = divider(signal.shape[0]) if q_oversample is None: - q_oversample = divider(n_coeff) + q_oversample, _ = divider(n_coeff) if window is None: window = np.exp(np.log(0.005) * np.linspace(-1, 1, nearest_odd(n_coeff)) ** 2) window = window / np.linalg.norm(window) diff --git a/tftb/tests/test_utils.py b/tftb/tests/test_utils.py index 47fc110..62801e5 100644 --- a/tftb/tests/test_utils.py +++ b/tftb/tests/test_utils.py @@ -28,11 +28,8 @@ def test_nextpow2(self): """Test the nextpow2 function.""" self.assertEqual(utils.nextpow2(2), 1) self.assertEqual(utils.nextpow2(17), 5) - import warnings - with warnings.catch_warnings(record=True) as catcher: + with self.assertRaises(ValueError): utils.nextpow2(-3) - self.assertEqual(len(catcher), 1) - self.assertTrue(catcher[-1].category, RuntimeWarning) def test_divider(self): """Test the divider function.""" diff --git a/tftb/utils.py b/tftb/utils.py index a9f4469..0045717 100644 --- a/tftb/utils.py +++ b/tftb/utils.py @@ -1,13 +1,8 @@ #! /usr/bin/env python -# -*- coding: utf-8 -*- -# vim:fenc=utf-8 -# -# Copyright © 2015 jaidev -# -# Distributed under terms of the MIT license. """Miscellaneous utilities.""" +import math import numpy as np @@ -61,9 +56,7 @@ def nextpow2(n): >>> print(nextpow2(x)) [ 0. 1. 2. 2. 3. 3. 3. 3.] """ - m_f = np.log2(n) - m_i = np.ceil(m_f) - return m_i + return math.ceil(math.log2(n)) def divider(N): @@ -78,20 +71,22 @@ def divider(N): :Example: >>> from __future__ import print_function >>> print(divider(256)) - (16.0, 16.0) + (16, 16) >>> print(divider(10)) - (2.0, 5.0) + (2, 5) >>> print(divider(101)) - (1.0, 101.0) + (1, 101) """ - n = np.floor(np.sqrt(N)) - while True: - old = n - m = np.ceil(N / float(n)) - n = np.floor(N / float(m)) - if n == old: + s = math.sqrt(N) + if s % 1 == 0: + s = int(s) + return s, s + s = math.ceil(s) + for i in range(s, 0, -1): + factor = N // i + if N % i == 0: break - return n, m + return i, factor def nearest_odd(N): @@ -109,18 +104,16 @@ def nearest_odd(N): >>> nearest_odd(3) 3.0 """ - if hasattr(N, "__iter__"): - N = np.array(N) + if isinstance(N, (list, np.ndarray)): y = np.floor(N) - y[np.remainder(y, 2) == 0] = np.ceil(N[np.remainder(y, 2) == 0]) y[np.remainder(y, 2) == 0] += 1 - return y + return y.astype(int) if N % 2 == 0: return N + 1 - elif np.floor(N) % 2 == 0: - return np.ceil(N) - elif np.floor(N) % 2 != 0: - return np.floor(N) + elif math.floor(N) % 2 == 0: + return math.ceil(N) + elif math.floor(N) % 2 != 0: + return math.floor(N) return N