% -------------------------------------------------------------------------------------------------
% 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<fmax))];
% Search the frequency of maximum power
[Pdmaxb, ibm] = max(Ppxx(ifband(1):ifband(2)));
Pbmax = 0.5*freq * Pdmaxb;
ibm = ifband(:,1)+ibm-1;
fPmaxb = Fp(ibm);
% ---------------------------------------------------------------------------------------------------------
% draw the marker of the power peak
fxband=[fPmaxb,fPmaxb];
ylimb=ylim+[0.1*ymima -0.1*ymima];
fyband = ones(size(fxband)); fyband(:,1)=ylimb(1); fyband(:,2)=ylimb(2);
plot(fxband',fyband', 'r-.');
% ---------------------------------------------------------------------------------------------------------
if(debug) fprintf('--- Press a key to write figure and results\n'); beep; pause; end
% Save on disk the figure
fig_dir='figure';
fig_path= fullfile(filepath_out, fig_dir);
if ~exist(fig_path,'dir') mkdir(filepath_out, fig_dir); end
fnamefig= fullfile(fig_path, cname);
%
set(gcf,'Color','white');
hTitle=14+ 8;
hXlabel=12+ 2;
hTlabel=hTitle+hXlabel;
bdwidth=5;
topbdwidth=60;
set(0,'Units','pixels')
scnsize = get(0,'ScreenSize');
fgwidth = scnsize(3)*.5-2*bdwidth;
fgheigth = scnsize(4)*.5;
pos = [bdwidth,...
scnsize(4)- fgheigth - (topbdwidth+bdwidth),...
fgwidth,...
fgheigth ];
set(gcf,'Position', pos);
cax = gca;
pos = get(cax,'Position');
hrTitle= hTitle/fgheigth;
hrXlabel= hXlabel/fgheigth;
% change the length of the main ordinate axis.
set(cax,'Position',[pos(1) pos(2)+hrXlabel pos(3) pos(4)-hrTitle]);
saveas(gcf,fnamefig, 'jpg')
% saveas(gcf,fnamefig, 'mfig')
% saveas(gcf,fnamefig, 'mmat')
% ---------------------------------------------------------------------------------------------------------
fpnamePar_out= fullfile(filepath_out, [cname '.mat']);
fprintf('Writing parameters and results on file: %s\n', fpnamePar_out);
save(fpnamePar_out, 'RRs', 'RRmean', 'RRstd', 'nc', 'crrf', 'rp1fn', 'rp2fn', ...
'fPmaxb', 'Pdmaxb');