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: 30]
        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:

- nt_bias_fft [c0,c1]=nt_bias_fft(x,freq,nfft) - covariance with and w/o filter bias
- nt_demean [y,mn]=nt_demean(x,w) - remove weighted mean over cols
- nt_dss0 [todss,pwr1,pwr2]=nt_dss0(c0,c1,keep1,keep2) - dss from covariance
- nt_greetings nt_greetings - display message the first time the toolbox is used
- nt_mmat y=nt_mmat(x,m) - matrix multiplication (with convolution)
- nt_pca [z,idx]=nt_pca(x,shifts,nkeep,threshold,w) - time-shift pca
- nt_smooth y=nt_smooth(x,T,nIterations,nodelayflag) - smooth by convolution with square window
- nt_spect_plot nt_spect_plot - plot power spectrum
- nt_tsr [y,idx,w]=nt_tsr(x,ref,shifts,wx,wref,keep,thresh) - time-shift regression (TSPCA)
- nt_wpwr [y,tweight]=nt_wpwr(x,w) - weighted power
- nt_zapline [y,yy]=nt_zapline(x,fline,nremove,p,plotflag) - remove power line artifact

This function is called by:

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

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