% ------------------------------------------------------------------------------------------------- % StftSpettriMedian.m: Spectral analysis of ecg signal after QRST complex canceling, % STFT is used. % This is an experimental routine to be tuned to improve atrial activity spectral estimation. % A more robust technique to reject artifactual spectra should be used before % computing the averaged spectrum. % % input parameter: % mainPath_in = main path for data files (generated by procedure "QrsAvgCancelTestClPcan.m") % casiDir = subdirectory with data files (learning-set, test-set-a, test-set-b) % cname = name of the case to process % ics = index of channel to process (usually 4) % freq = sampling frequency % input files (generated by the procedure "QrsAvgCancelTestClPcan.m"): % signal file (.dat), usually 4 channels the last obtained cancelin QRST from ecg signals. % parameter file (.mat), containing: 'RRs', 'RRmean', 'RRstd', 'nc', 'crrf', 'rp1fn', 'rp2fn' % output files: % (directory "SptRes") parameter file (.mat) updated with spectral measures: 'fPmaxb', 'Pdmaxb' % (directory "SptRes\figure") figure file (.jpg) representing the spectrum of the estimate atrial signal % % Copyright (C) 2004 Maurizio Varanini, Clinical Physiology Institute, CNR, Pisa, Italy % % This program is free software; you can redistribute it and/or modify it under the terms % of the GNU General Public License as published by the Free Software Foundation; either % version 2 of the License, or (at your option) any later version. % % This program is distributed "as is" and "as available" in the hope that it will be useful, % but WITHOUT ANY WARRANTY of any kind; without even the implied warranty of MERCHANTABILITY % or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. % % You should have received a copy of the GNU General Public License along with this program; % if not, write to the Free Software Foundation, Inc., % 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. % % For any comment or bug report, please send e-mail to: maurizio.varanini@ifc.cnr.it % ------------------------------------------------------------------------------------------------- function StftSpettriMedian(mainPath_in, casiDir, cname, ics, freq) %nargin=0 echo off % clc close all debug=0; graf=0; SWfilt=0; datiModExt=''; if nargin<1 % clear all debug=0; graf=1; SWfilt=0; mainPath_in='..\..\DatiQcD'; datiModExt=''; freq=1024; % casiDir='learning-set'; casiDir='test-set-a'; % cname='t06'; cname='a10'; % cname='a29'; % cname=minput('Name of case to process','%s',cname); end [progpath, progname] = fileparts(which(mfilename)); fprintf('\n-------------------------------------------------------------------\n'); fprintf('Name of case to process: %s\n', cname); mainPath_in= fullfile(mainPath_in, casiDir,''); filepath=mainPath_in; extension=[datiModExt, '.dat']; ns=4; fpnamePar_in=fullfile(filepath, [cname '.mat']); fprintf('Lettura parametri dal file: %s\n', fpnamePar_in); load(fpnamePar_in); % RRs RRmean RRstd nc crrf rp1fn rp2fn % fpname= strcat(filepath, cname, extension); fpname = fullfile(filepath, [cname extension]); fprintf('Lettura del file: %s\n', fpname); fid = fopen(fpname,'r'); [datap, count] = fread(fid,[ns,inf],'double'); fclose(fid); % case type: 1->N Not terminated, 2->T Terminated tipoC=cname(1); %------------------------------------------------------------- if nargin<3 % index of signal to process ics=4; % ics=minput('indice segnale da elaborare','%d',ics); end %------------------------------------------------------------- dirOut='SptRes'; if debug dirOut=[dirOut 'D']; end fprintf('case type= %s, signal index= %d\n', tipoC, ics); fprintf('Directory for writing results = %s\n', dirOut); filepath_out= fullfile(mainPath_in, dirOut); if ~exist(filepath_out,'dir') mkdir(mainPath_in, dirOut); end if(freq > 256) % decimation to a frequency close to 256Hz stpd=round(freq/256); data=datap(:,1:stpd:end); freq=freq/stpd; else data=datap(:,1:end-freq); end % registration length [ns ndt]=size(data); lendts=ndt/freq; fprintf('Sampling frequency = %g\n', freq); fprintf('Total recording length = %g (samples), = %g (seconds)\n', ndt, lendts); ysgn=data(ics,:)'; if(graf) % plot of signal PlotSgnMrkN(ysgn, [], freq, [cname, ' signal-', num2str(ics)]); end iniint=1; finint=ndt; % finint=length(ysgn)- 1*freq; % iniint=finint- 10*freq; nd= finint - iniint; lends=nd/freq; fprintf('begin - end of signal interval in samples= %g - %g\n', iniint, finint); fprintf('interval length= %g (samples), = %g (seconds)\n', nd, lends); x=data(ics,iniint:finint); % x=detrend(x); % figure; % plot(x,'g'); %hold on if(SWfilt) fmind=3; fmaxd=10; orderFir=fix(3*freq/fmind); Wn = [fmind fmaxd]/(freq/2); % The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0 % corresponding to half the sample rate B = FIR1(orderFir,Wn); % default -> Hamming fprintf('number of filter coeff.= %d', length(B)); x = filter(B,1,x); % plot(x,'b'); hold off; % title([cname ' - signal and filtered signal']); if(graf) PlotSgnMrkN(x, [], freq, [cname, ' filtered signal-', num2str(ics)]); end end % --------------------------------------------------------------------------------------------------------- % --------------------------------------------------------------------------------------------------------- % dfmin=0.25; % dfmin=0.15; dfmin=0.2; % frequency resolution WLEN=fix(1/dfmin*freq); NFFT=2.^(fix(log2(WLEN))+2); if (WLEN>nd) WLEN=nd; end OVERLAP=fix(0.95 * WLEN); risFreq= freq/WLEN; % [STFT,Fp,T] = specgram(x,2048,freq,hann(512),500); [STFT,Fp,T] = specgram(x,NFFT,freq,gausswin(WLEN,3.5),OVERLAP); %[STFT,Fp,T] = specgram(x,NFFT,freq,hamming(WLEN),OVERLAP); if(graf) figure; imagesc(T,Fp,20*log10(abs(STFT))); axis xy, colormap(JET); end % each column of STFT contains an estimate of the short-term, time-localized frequency content of the signal ew=gausswin(WLEN,3.5)'*gausswin(WLEN,3.5)/NFFT; % ew=hamming(WLEN)'*hamming(WLEN)/NFFT; Smagn= abs(STFT); Spow= (2/NFFT) *(Smagn .* Smagn)/ew; Spow(1,:)= Spow(1,:)/2; if (rem(NFFT,2)==0) Spow(end,:)= Spow(end,:)/2; end Spow= Spow/freq; % --------------------------------------------------------------------------------------------------------- Ppxx=mean(Spow,2); % Check total power Pptot = (0.5*freq/length(Ppxx)) * sum(Ppxx); Pxtot= x*x'/length(x); fprintf('Total power Pxtot= %g, Pptot= %g\n', Pxtot, Pptot); % Attempt to build a robust PSD estimate rejecting "wild" values for i=1:size(Spow,1) Ppxx(i,1)=meansc(Spow(i,:),10,10); end if(graf) figure; psdplot(Ppxx,Fp,'Hz','db',[cname 'PSD Estimate']); ylim = get(gca,'Ylim'); if(debug) fprintf('--- Press a key to continue\n'); pause; end end % --------------------------------------------------------------------------------------------------------- fmin=3; fmax=10; if(graf) figure psdplot(Ppxx,Fp,'Hz','db',[cname ' - PSD Estimate']); xlim=[fmin fmax]; set(gca,'Xlim',xlim); set(gca,'Xscale','log'); ylim = get(gca,'Ylim'); hold on end % --------------------------------------------------------------------------------------------------------- figure psdplot(Ppxx,Fp,'Hz','db',[cname ' - PSD Estimate']); hold on xlim=[fmin fmax]; set(gca,'Xlim',xlim); set(gca,'Xscale','linear'); ylim = get(gca,'Ylim'); ymima=ylim(2)-ylim(1); % --------------------------------------------------------------------------------------------------------- % draw the marker of decision threshold frequency fThreshold=5.5; fxband=[fThreshold,fThreshold]; ylimb=ylim+[0.0*ymima -0.0*ymima]; fyband = ones(size(fxband)); fyband(:,1)=ylimb(1); fyband(:,2)=ylimb(2); plot(fxband',fyband', 'g--'); % --------------------------------------------------------------------------------------------------------- ifband= [min(find(Fp>fmin)), max(find(Fp