function varargout=svdslep(N,NW,K,method,imp) % [E,V,SE,ngro]=SVDSLEP(N,NW,K,method,imp) % % Explicit diagonalization of the Slepian concentration operators % in one linear dimension. % % INPUT: % % N The length of the sequence [default: 2^8] % NW Half the Shannon number (time-bandwidth product) [default: 10] % K The requested number of eigenvectors [default: 2NW-1] % method 1 Using EIGS on the Hermitian form [default] % 2 Using EIG on the Hermitian form % 3 Using SVDS on the projection operator % 4 Using SVD on the projection operator % 5 Using the power method for the largest eigenvalue % imp 1 Uses the implicit operator method [default] % 0 Uses the sparse explicit operator method [slow] % % OUTPUT: % % E The eigenfunctions of the concentration problem % V The eigenvalues of the concentration problem % SE The power spectrum of the eigenfunctions % ngro The "growth" factor used % % EXAMPLE: % % svdslep('demo1') % % Last modified by fjsimons-at-alum.mit.edu, 06/21/2010 % See some notes in RB VIII p 90 % Default values defval('N',2^10) if ~isstr(N) defval('NW',10) defval('K',2*NW-1) defval('method',1) defval('ngro',256); defval('imp',1); % This one for "excessive" verification defval('xver',0) % Make a larger domain Nd=N*ngro; if ngro>32 || N>256 % Force to the implicit method or your machine will die imp=1; disp('Method reset to IMPLICIT') end if imp==1 disp(sprintf('Implicit method embedded %ix\n',ngro)) nonz=(Nd-N)/2+[1:N]; P =@(x) proj(x,nonz); Q =@(x) fft(x); Qi=@(x) ifft(x); nonz=[1:ngro*NW+1 Nd-ngro*NW+1:Nd]; L =@(y) proj(y,nonz); H =@(x) P(Qi(L(Q(P(x))))); % Acknowledge that H is complex (though it is symmetric) OPTS.isreal=false; OPTS.disp=0; % Remember to specify the output size [E,V]=eigs(H,Nd,K,'LR',OPTS); [V,i]=sort(diag(V),'descend'); E=E(:,i); V=V(1:K); E=E(:,1:K); elseif imp==0 disp(sprintf('Explicit, sparse method embedded %ix\n',ngro)) % The spatial projection operator P=diag([zeros(1,(Nd-N)/2) ones(1,N) zeros(1,(Nd-N)/2)]); % The Fourier transform operator Q=dftmtx(Nd)/sqrt(Nd); % The bandlimiting operator L=diag([ones(1,ngro*NW+1) zeros(1,Nd-2*ngro*NW-1) ones(1,ngro*NW)]); % Make them all sparse for speed - could start from it altogether P=sparse(P); Q=sparse(Q); L=sparse(L); if xver==1 % Check that Q is a unitary matrix difer(Q'*Q-diag(diag(Q'*Q))) difer(Q'*Q-eye(size(Q))) % Check that H is a Hermitian matrix difer(H'-H) % Check that H is a positive definite matrix [R,p]=chol(H); isposdef=p==0; difer(isposdef-1) % Check that Q does what I think it does f=rand(Nd,1); % See the normalization of FFT (1) and IFFT (1/N) difer(fft(f)/sqrt(Nd)-Q*f); end % The operator whose singular functions we want A=L*Q*P; % The operator whose eigenfunctions we want H=P*Q'*A; % The singular values of A are the eigenvalues of AA'. The singular % values of the concentration/projection operator are the eigenvalues of % the "squared" projection operator. % The eigenvector decomposition switch method case 1 disp('Using EIGS') OPTS.disp=0; [E,V]=eigs(H,K,'LR',OPTS); [V,i]=sort(diag(V),'descend'); E=E(:,i); V=V(1:K); E=E(:,1:K); case 2 disp('Using EIG') [E,V]=eig(H,'nobalance'); [V,i]=sort(diag(V),'descend'); E=E(:,i); V=V(1:K); E=E(:,1:K); case 3 disp('Using SVDS') [U,V,E]=svds(A,K); [V,i]=sort(diag(V).^2,'descend'); E=E(:,i); V=V(1:K); E=E(:,1:K); case 4 disp('Using SVD') [U,V,E]=svd(A); [V,i]=sort(diag(V).^2,'descend'); E=E(:,i); V=V(1:K); E=E(:,1:K); case 5 disp('Using the power method for the first one only') Vi=1; Vj=2; E=rand(ngro*N,1); E=E/sqrt(E'*E); while abs(Vj-Vi)>1e-6 Vi=Vj; Vj=E; E=H*E; Vj=real(E'*Vj); E=real(E/sqrt(E'*E)); plot(E); pause(0.1) end V=Vj; end end % Define some kind of tolerance level tol=100*eps; % Make them real as we know they should be if any(imag(V)>tol) error('Complex eigenvalues'); else V=real(V); % Check imaginary part of the "good" eigenfunctions disp(sprintf('mean(abs(imag(E))) = %8.3e out of %8.3e\n',... mean(mean(abs(imag(E(:,V>tol))))),... mean(mean(abs(E(:,V>tol)))))) % Take out only the central part E=E((Nd-N)/2+1:(Nd-N)/2+N,:); % Note that they were normalized in the complex plane E=real(E); E=E./repmat(diag(sqrt(E'*E))',size(E,1),1); % Cannot do abs(E).*sign(E) because the sign is bad for complex % But either real or imag seems to work... they more or less scale! end if xver==1 % Check the orthogonality for the "good" eigenfunctions ortho=E(:,V>tol)'*E(:,V>tol); difer(ortho-diag(diag(ortho))) difer(diag(ortho)-1) end if nargout>2 % Get the power spectrum SE=abs(fft(E,8*length(E))).^2; else SE=NaN; end % Fix the first bit to be an upswing for index=1:size(E,2) E(:,index)=E(:,index)*sign(E(2,index)); end % Output varns={E,V,SE,ngro}; varargout=varns(1:nargout); elseif strcmp(N,'demo1') % So you can say svdslep('demo1',1) defval('NW',[]) imp=NW; defval('imp',1) N=2^(3+ceil(rand*5)); NW=ceil(rand*5); N=2^8; NW=3; method=1; disp(sprintf('\nN = %i , NW = %i\n',N,NW)) [E,V,SE,ngro]=svdslep(N,NW,max(3,2*NW),method,imp); % Compare with Matlab's built-in functions [E2,V2]=dpss(N,NW,max(3,2*NW)); clf [ah,ha]=krijetem(subnum(2,3)); % Make the frequency axis for the positive frequencies only [fax,selekt]=fftaxis1d(E(:,1),length(SE),N-1); SE=SE(selekt,:); % Another, identical, more recently revised, way % [K,kx]=knum2([2 N],[2 N]-1); % fx=-fliplr(indeks(kx,selekt)/2/pi); yls=[-100 2]; for ind=1:3 SE2=indeks(abs(fft(E2(:,ind),8*length(E2)).^2),selekt); axes(ah(ind)) plot(E(:,ind),'b','linew',2); hold on plot(E2(:,ind),'y') disp(sprintf('Eigenfunction %i: rms %8.3f%s',ind,... rms(E(:,ind)-E2(:,ind))/rms(E(:,ind))*100,'%')) yl(ind)=ylabel('Slepian eigenfunction'); set(ah(ind),'xtick',[1 N],'xlim',[1 N]) xlabel(sprintf('%s = %12.9f (EIGS)','\lambda',V(ind))) axes(ha(2*ind)) plot(N*fax,decibel(SE(:,ind)),'b','linew',2); hold on plot(N*fax,decibel(SE2),'y'); plot(N*[NW/N NW/N],yls,'k-'); plot(-N*[NW/N NW/N],yls,'k-'); hold off %xl(3+ind)=xlabel('frequency (sample)^{-1}'); xl(3+ind)=xlabel('frequency x sample size'); yl(3+ind)=ylabel('power spectral density (dB)'); title(sprintf('%s = %12.9f (DPSS)','\lambda',V2(ind))) % Show the entire range too=8*1.95; exes=[0 1/too]; exes=[-0.1 20]; set(ha(2*ind),'xtick',round(exes*10)/10,'xlim',minmax(exes),... 'xgrid','off','ylim',yls) end % Cosmetics delete(yl([2 3 5 6])) longticks(ah) tt=supertit(ah(1:3),... sprintf('SVDSLEP vs DPSS on a %ix domain',ngro),12); movev(tt,0.15) fig2print(gcf,'landscape') % Is the ratio of the two things symmetric? % for ind=1:size(E,2) % a=E(:,ind)./E2(:,ind); % rmsym(ind)=rms([a(1:length(a)/2)-flipud(a(length(a)/2+1:end))])./... % rms(a)*100; % end % rmsym figdisp([],'demo1') end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Pv=proj(v,p) % Projects the vector v on the indices p Pv=zeros(size(v)); Pv(p)=v(p);