Home > NoiseTools > nt_zapline.m

nt_zapline

PURPOSE ^

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

SYNOPSIS ^

function [y,yy]=nt_zapline(x,fline,nremove,p,plotflag)

DESCRIPTION ^

[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:1024]
    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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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:1024]
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=1024; 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=size(x,1); 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('original'); 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(y/sqrt(mean(x(:).^2)),p.nfft,[],[],1/fline);
0067     semilogy(f,abs(pxx)/divisor);
0068     hold on
0069     [pxx,f]=nt_spect_plot((x-y)/sqrt(mean(x(:).^2)),p.nfft,[],[],1/fline);
0070     semilogy(f,abs(pxx)/divisor);
0071     legend('clean', 'removed'); legend boxoff
0072     set(gca,'ygrid','on','xgrid','on');
0073     set(gca,'yticklabel',[]); ylabel([]);
0074     xlabel('frequency (relative to line)');
0075     yl2=get(gca,'ylim');
0076     hh=get(gca,'children');
0077     set(hh(1),'color',[1 .5 .5]); set(hh(2), 'color', [ 0 .7 0]); 
0078     set(hh(2),'linewidth', 2);
0079       
0080        yl(1)=min(yl1(1),yl2(1)); yl(2)=max(yl1(2),yl2(2));
0081     subplot 121; ylim(yl); subplot 122; ylim(yl);
0082     drawnow;
0083     return
0084 end
0085 if ~nargout
0086     clear y yy
0087     return
0088 end
0089 
0090 xx=nt_smooth(x,1/fline,p.niterations); % cancels line_frequency and harmonics, light lowpass
0091 if isempty(p.nkeep); p.nkeep=size(x,2); end
0092 xxxx=nt_pca(x-xx,[],p.nkeep); % reduce dimensionality to avoid overfitting
0093 
0094 % DSS to isolate line components from residual:
0095 nHarmonics=floor((1/2)/fline);
0096 [c0,c1]=nt_bias_fft(xxxx,fline*(1:nHarmonics), p.nfft);
0097 
0098 [todss,pwr0,pwr1]=nt_dss0(c0,c1);
0099 assert(size(todss,2)>0, '!'); 
0100 if ~nargout;
0101     figure(p.fig1); clf;
0102     plot(pwr1./pwr0, '.-'); xlabel('component'); ylabel('score'); title('DSS to enhance line frequencies');
0103 end
0104 xxxx=nt_mmat(xxxx,todss(:,1:nremove)); % line-dominated components
0105 xxx=nt_tsr(x-xx,xxxx); % project them out
0106 clear xxxx
0107 
0108 % reconstruct clean signal
0109 y=xx+xxx; clear xx xxx
0110 yy=x-y;
0111 
0112 % test code
0113 if 0
0114     sr=400;
0115     nsamples=100000; nchans=100;
0116     signal=randn(nsamples,nchans);
0117     artifact=sin((1:nsamples)'/sr*2*pi*50);
0118     artifact=max(artifact,0).^3; % introduce harmonics
0119     artifact=3*nt_demean(artifact*randn(1,nchans));
0120     disp(nt_wpwr(artifact)/nt_wpwr(signal+artifact));
0121     nt_zapline(signal+artifact,50/sr);
0122 end
0123

Generated on Thu 07-May-2020 08:42:54 by m2html © 2005