0001 function [y,yy]=nt_zapline(x,fline,nremove,p,plotflag)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
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
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);
0091 if isempty(p.nkeep); p.nkeep=size(x,2); end
0092 xxxx=nt_pca(x-xx,[],p.nkeep);
0093
0094
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));
0105 xxx=nt_tsr(x-xx,xxxx);
0106 clear xxxx
0107
0108
0109 y=xx+xxx; clear xx xxx
0110 yy=x-y;
0111
0112
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;
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