nt_narrowband_scan

A=nt_narrowband_scan(x,freqs,sr,Q,plotflag) - scan for narrowband components using DSS

  A: cell array of DSS matrices

  x: data (time*channels or time*channels*trials)
  freqs: Hz, array of bias frequencies
  sr: Hz, sampling rate
  Q: quality factors of scanning filter (default: [8 4])
  plotflag: if true plot (default: 1)

 If no output arguments, plots spectra of first DSS components for each
 bias

0001 function A=nt_narrowband_scan(x,freqs,sr,Q,plotflag)
0002 %A=nt_narrowband_scan(x,freqs,sr,Q,plotflag) - scan for narrowband components using DSS
0003 %
0004 %  A: cell array of DSS matrices
0005 %
0006 %  x: data (time*channels or time*channels*trials)
0007 %  freqs: Hz, array of bias frequencies
0008 %  sr: Hz, sampling rate
0009 %  Q: quality factors of scanning filter (default: [8 4])
0010 %  plotflag: if true plot (default: 1)
0011 %
0012 % If no output arguments, plots spectra of first DSS components for each
0013 % bias
0014 %
0015 nt_greetings;
0016 
0017 if nargin<5; plotflag=[]; end
0018 if nargin<4||isempty(Q); Q=[8 4]; end
0019 if nargin<3; error('!'); end
0020 
0021 freqs=freqs(:);
0022 
0023 
0024 current_figure=get(0,'CurrentFigure');
0025 A={};
0026 for iFreqs=1:numel(freqs)
0027     freq=freqs(iFreqs);
0028     [b,a]=nt_filter_peak(freq/(sr/2),Q(1));
0029     [c0,c1]=nt_bias_filter(x,b,a);
0030     if numel(Q)==2
0031         [b,a]=nt_filter_peak(freq/(sr/2),Q(2));
0032         [~,c0]=nt_bias_filter(x,b,a);
0033     end
0034     [todss,pwr0,pwr1]=nt_dss0(c0,c1);
0035     
0036     if 1; figure(100); clf; plot(pwr1./pwr0,'.-'); title([num2str(freq), 'Hz bias']); ylabel('score');  drawnow;end
0037     
0038     A{iFreqs}=todss;
0039 end
0040 
0041 
0042 if ~nargout || ~isempty(plotflag)
0043     if isempty(current_figure);
0044         figure; 
0045     else
0046         figure(current_figure);
0047     end
0048     AA=zeros(size(x,2),numel(freqs));
0049     for iFreqs=1:numel(freqs)
0050         AA(:,iFreqs)=A{iFreqs}(:,1);
0051     end
0052     x=nt_mmat(x,AA);
0053     nfft=2.^nextpow2(size(x,1)+1)/2;
0054     MAX_NFFT=1024;  
0055     nfft=min(nfft,MAX_NFFT);
0056     nt_spect_plot2(nt_normcol(x),nfft,[],[],sr);
0057     K=round(numel(freqs)/6);
0058     set(gca,'ytick',1:K:numel(freqs), 'yticklabel',num2str(freqs(1:K:end), '%.3g')); 
0059     ylabel('bias, Hz'); xlabel('frequency, Hz');
0060     set(gca,'xgrid','on','xminortick','on');
0061     drawnow;
0062 end
0063 
0064 if ~nargout;     clear A; end