nt_zapline

[y,yy]=nt_zapline(x,fline,nremove,p,plotflag) - remove power line artifact

  y: denoised data
  yy: artifact

  x: data
  fline: line frequency (normalized to sr)
  nremove: number of components to remove [default: 1]
  p: additional parameters:
    p.nfft: size of FFT [default:512]
    p.nkeep: number of components to keep in DSS [default: all]
    p.niterations: number of iterations for smoothing filter
    p.fig1: figure to use for DSS score [default: 100]
    p.fig2: figure to use for results [default: 101]
  plotflag: plot

Examples:
  nt_zapline(x,60/1000) 
    apply to x, assuming line frequency=60Hz and sampling rate=1000Hz, plot results
  nt_zapline(x,60/1000,4)
    same, removing 4 line-dominated components 
  p=[];p.nkeep=30; nt_zapline(x,60/1000,4,p);
    same, truncating PCs beyond the 30th to avoid overfitting
  [y,yy]=nt_zapline(x,60/1000)
    return cleaned data in y, noise in yy, don't plot

0001 function [y,yy]=nt_zapline(x,fline,nremove,p,plotflag)
0002 %[y,yy]=nt_zapline(x,fline,nremove,p,plotflag) - remove power line artifact
0003 %
0004 %  y: denoised data
0005 %  yy: artifact
0006 %
0007 %  x: data
0008 %  fline: line frequency (normalized to sr)
0009 %  nremove: number of components to remove [default: 1]
0010 %  p: additional parameters:
0011 %    p.nfft: size of FFT [default:512]
0012 %    p.nkeep: number of components to keep in DSS [default: all]
0013 %    p.niterations: number of iterations for smoothing filter
0014 %    p.fig1: figure to use for DSS score [default: 100]
0015 %    p.fig2: figure to use for results [default: 101]
0016 %  plotflag: plot
0017 %
0018 %Examples:
0019 %  nt_zapline(x,60/1000)
0020 %    apply to x, assuming line frequency=60Hz and sampling rate=1000Hz, plot results
0021 %  nt_zapline(x,60/1000,4)
0022 %    same, removing 4 line-dominated components
0023 %  p=[];p.nkeep=30; nt_zapline(x,60/1000,4,p);
0024 %    same, truncating PCs beyond the 30th to avoid overfitting
0025 %  [y,yy]=nt_zapline(x,60/1000)
0026 %    return cleaned data in y, noise in yy, don't plot
0027 %
0028 
0029 % NoiseTools
0030 try, nt_greetings; catch, disp('You must download NoiseNools from http://audition.ens.fr/adc/NoiseTools/'); return; end
0031 
0032 assert(nargin>=2, '!'); 
0033 if nargin<3||isempty(nremove); nremove=1; end
0034 if nargin<4; p=[]; end
0035 if ~isfield(p,'nfft'); p.nfft=512; end
0036 if ~isfield(p,'nkeep'); p.nkeep=[]; end
0037 if ~isfield(p,'niterations'); p.niterations=1; end
0038 if ~isfield(p,'fig1'); p.fig1=100; end
0039 if ~isfield(p, 'fig2'); p.fig2=101; end
0040 if nargin<5||isempty(plotflag); plotflag=0; end
0041 
0042 if isempty(x); error('!'); end
0043 if nremove>=size(x,1); error('!'); end
0044 if fline>1/2; error('fline should be less than Nyquist'); end
0045 if size(x,1)<p.nfft; warning(['reducing nfft to ',num2str(size(x,1))]); p.nfft=2*floor(size(x,1)/2); end
0046 
0047 if ~nargout || plotflag
0048     % print result and display spectra
0049     y=nt_zapline(x,fline,nremove,p);
0050     disp('proportion of non-DC power removed:');
0051     disp(nt_wpwr(x-y)/nt_wpwr(nt_demean(x)));
0052     
0053     figure(p.fig2); clf;    
0054     subplot 121
0055     [pxx,f]=nt_spect_plot(x/sqrt(mean(x(:).^2)),p.nfft,[],[],1/fline);
0056     divisor=sum(pxx);
0057     semilogy(f,abs(pxx)/divisor);
0058     legend('raw'); legend boxoff
0059     set(gca,'ygrid','on','xgrid','on');
0060     xlabel('frequency (relative to line)');
0061     ylabel('relative power');
0062     yl1=get(gca,'ylim');
0063     hh=get(gca,'children');
0064     set(hh(1),'color','k')
0065     subplot 122
0066     [pxx,f]=nt_spect_plot((x-y)/sqrt(mean(x(:).^2)),p.nfft,[],[],1/fline);
0067     semilogy(f,abs(pxx)/divisor);
0068     hold on
0069     [pxx,f]=nt_spect_plot(y/sqrt(mean(x(:).^2)),p.nfft,[],[],1/fline);
0070     semilogy(f,abs(pxx)/divisor);
0071     hold on
0072     legend('noise', 'clean'); legend boxoff
0073     set(gca,'ygrid','on','xgrid','on');
0074     set(gca,'yticklabel',[]); ylabel([]);
0075     xlabel('frequency (relative to line)');
0076     yl2=get(gca,'ylim');
0077     hh=get(gca,'children');
0078     set(hh(2),'color',[1 .5 .5]); 
0079     set(hh(1), 'linewidth', 0.25);
0080     set(hh(1), 'color', [ 0 .7 0]); 
0081     set(hh(1),'linewidth', 2);
0082       
0083        yl(1)=min(yl1(1),yl2(1)); yl(2)=max(yl1(2),yl2(2));
0084     subplot 121; ylim(yl); subplot 122; ylim(yl);
0085     drawnow;
0086     return
0087 end
0088 if ~nargout
0089     clear y yy
0090     return
0091 end
0092 
0093 xx=nt_smooth(x,1/fline,p.niterations); % cancels line_frequency and harmonics, light lowpass
0094 if isempty(p.nkeep); p.nkeep=size(x,2); end
0095 xxxx=nt_pca(x-xx,[],p.nkeep); % reduce dimensionality to avoid overfitting
0096 
0097 % DSS to isolate line components from residual:
0098 nHarmonics=floor((1/2)/fline);
0099 [c0,c1]=nt_bias_fft(xxxx,fline*(1:nHarmonics), p.nfft);
0100 
0101 [todss,pwr0,pwr1]=nt_dss0(c0,c1);
0102 assert(size(todss,2)>0, '!'); 
0103 if ~nargout;
0104     figure(p.fig1); clf;
0105     plot(pwr1./pwr0, '.-'); xlabel('component'); ylabel('score'); title('DSS to enhance line frequencies');
0106 end
0107 xxxx=nt_mmat(xxxx,todss(:,1:nremove)); % line-dominated components
0108 xxx=nt_tsr(x-xx,xxxx); % project them out
0109 clear xxxx
0110 
0111 % reconstruct clean signal
0112 y=xx+xxx; clear xx xxx
0113 yy=x-y;
0114 
0115 % test code
0116 if 0
0117     sr=400;
0118     nsamples=100000; nchans=100;
0119     signal=randn(nsamples,nchans);
0120     artifact=sin((1:nsamples)'/sr*2*pi*50);
0121     artifact=max(artifact,0).^3; % introduce harmonics
0122     artifact=3*nt_demean(artifact*randn(1,nchans));
0123     disp(nt_wpwr(artifact)/nt_wpwr(signal+artifact));
0124     nt_zapline(signal+artifact,50/sr);
0125 end
0126