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.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.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 isempty(x); error('!'); end
0041 if nremove>=size(x,1); error('!'); end
0042 if fline>1/2; error('fline should be less than Nyquist'); end
0043 if size(x,1)<p.nfft; warning(['reducing nfft to ',num2str(size(x,1))]); p.nfft=size(x,1); end
0044 
0045 if ~nargout || plotflag
0046     % print result and display spectra
0047     y=nt_zapline(x,fline,nremove,p);
0048     disp('proportion of non-DC power removed:');
0049     disp(nt_wpwr(x-y)/nt_wpwr(nt_demean(x)));
0050     
0051     figure(p.fig2); clf;    
0052     subplot 121
0053     [pxx,f]=nt_spect_plot(x/sqrt(mean(x(:).^2)),p.nfft,[],[],1/fline);
0054     divisor=sum(pxx);
0055     semilogy(f,abs(pxx)/divisor);
0056     legend('original'); legend boxoff
0057     set(gca,'ygrid','on','xgrid','on');
0058     xlabel('frequency (relative to line)');
0059     ylabel('relative power');
0060     yl1=get(gca,'ylim');
0061     hh=get(gca,'children');
0062     set(hh(1),'color','k')
0063     subplot 122
0064     [pxx,f]=nt_spect_plot(y/sqrt(mean(x(:).^2)),p.nfft,[],[],1/fline);
0065     semilogy(f,abs(pxx)/divisor);
0066     hold on
0067     [pxx,f]=nt_spect_plot((x-y)/sqrt(mean(x(:).^2)),p.nfft,[],[],1/fline);
0068     semilogy(f,abs(pxx)/divisor);
0069     legend('clean', 'removed'); legend boxoff
0070     set(gca,'ygrid','on','xgrid','on');
0071     set(gca,'yticklabel',[]); ylabel([]);
0072     xlabel('frequency (relative to line)');
0073     yl2=get(gca,'ylim');
0074     hh=get(gca,'children');
0075     set(hh(1),'color',[1 .5 .5]); set(hh(2), 'color', [ 0 .7 0]); 
0076     set(hh(2),'linewidth', 2);
0077       
0078        yl(1)=min(yl1(1),yl2(1)); yl(2)=max(yl1(2),yl2(2));
0079     subplot 121; ylim(yl); subplot 122; ylim(yl);
0080     drawnow;
0081     return
0082 end
0083 if ~nargout
0084     clear y yy
0085     return
0086 end
0087 
0088 xx=nt_smooth(x,1/fline); % cancels line_frequency and harmonics, light lowpass
0089 if isempty(p.nkeep); p.nkeep=size(x,2); end
0090 xxxx=nt_pca(x-xx,[],p.nkeep); % reduce dimensionality to avoid overfitting
0091 
0092 % DSS to isolate line components from residual:
0093 nHarmonics=floor((1/2)/fline);
0094 [c0,c1]=nt_bias_fft(xxxx,fline*(1:nHarmonics), p.nfft);
0095 
0096 [todss,pwr0,pwr1]=nt_dss0(c0,c1);
0097 assert(size(todss,2)>0, '!'); 
0098 if ~nargout;
0099     figure(p.fig1); clf;
0100     plot(pwr1./pwr0, '.-'); xlabel('component'); ylabel('score'); title('DSS to enhance line frequencies');
0101 end
0102 xxxx=nt_mmat(xxxx,todss(:,1:nremove)); % line-dominated components
0103 xxx=nt_tsr(x-xx,xxxx); % project them out
0104 clear xxxx
0105 
0106 % reconstruct clean signal
0107 y=xx+xxx; clear xx xxx
0108 yy=x-y;
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

Generated on Tue 18-Feb-2020 11:23:12 by m2html © 2005