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 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
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);
0089 if isempty(p.nkeep); p.nkeep=size(x,2); end
0090 xxxx=nt_pca(x-xx,[],p.nkeep);
0091
0092
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));
0103 xxx=nt_tsr(x-xx,xxxx);
0104 clear xxxx
0105
0106
0107 y=xx+xxx; clear xx xxx
0108 yy=x-y;
0109
0110
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;
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