[R] Conversion of Matlab code to an R code
John C Frain
frainj at gmail.com
Mon Mar 23 21:28:39 CET 2015
You do not say why you want to convert. I your reason is that you do not
have access to Matlab you might try to run your program in Octave. I think
that there are some syntax errors in your Matlab code. Most of these, but
not all, I think, may have been caused, (as you have already been told) by
the use of HTML in your email.
John C Frain
3 Aranleigh Park
Rathfarnham
Dublin 14
Ireland
www.tcd.ie/Economics/staff/frainj/home.html
mailto:frainj at tcd.ie
mailto:frainj at gmail.com
On 23 March 2015 at 15:10, Abhinaba Roy <abhinabaroy09 at gmail.com> wrote:
> Hi,
>
> Can a Matlab code be converted to R code?
>
> I am finding it difficult to do so.
>
> Could you please help me out with it.
>
> Your help will be highly appreciated.
>
> Here comes the Matlab code
> ##############
> if ~isvector(ecg)
> error('ecg must be a row or column vector');
> end
>
>
> if nargin < 3
> gr = 1; % on default the function always plots
> end
> ecg = ecg(:); % vectorize
>
> %% Initialize
> qrs_c =[]; %amplitude of R
> qrs_i =[]; %index
> SIG_LEV = 0;
> nois_c =[];
> nois_i =[];
> delay = 0;
> skip = 0; % becomes one when a T wave is detected
> not_nois = 0; % it is not noise when not_nois = 1
> selected_RR =[]; % Selected RR intervals
> m_selected_RR = 0;
> mean_RR = 0;
> qrs_i_raw =[];
> qrs_amp_raw=[];
> ser_back = 0;
> test_m = 0;
> SIGL_buf = [];
> NOISL_buf = [];
> THRS_buf = [];
> SIGL_buf1 = [];
> NOISL_buf1 = [];
> THRS_buf1 = [];
>
>
> %% Plot differently based on filtering settings
> if gr
> if fs == 200
> figure, ax(1)=subplot(321);plot(ecg);axis tight;title('Raw ECG Signal');
> else
> figure, ax(1)=subplot(3,2,[1 2]);plot(ecg);axis tight;title('Raw ECG
> Signal');
> end
> end
> %% Noise cancelation(Filtering) % Filters (Filter in between 5-15 Hz)
> if fs == 200
> %% Low Pass Filter H(z) = ((1 - z^(-6))^2)/(1 - z^(-1))^2
> b = [1 0 0 0 0 0 -2 0 0 0 0 0 1];
> a = [1 -2 1];
> h_l = filter(b,a,[1 zeros(1,12)]);
> ecg_l = conv (ecg ,h_l);
> ecg_l = ecg_l/ max( abs(ecg_l));
> delay = 6; %based on the paper
> if gr
> ax(2)=subplot(322);plot(ecg_l);axis tight;title('Low pass filtered');
> end
> %% High Pass filter H(z) = (-1+32z^(-16)+z^(-32))/(1+z^(-1))
> b = [-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 32 -32 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 1];
> a = [1 -1];
> h_h = filter(b,a,[1 zeros(1,32)]);
> ecg_h = conv (ecg_l ,h_h);
> ecg_h = ecg_h/ max( abs(ecg_h));
> delay = delay + 16; % 16 samples for highpass filtering
> if gr
> ax(3)=subplot(323);plot(ecg_h);axis tight;title('High Pass Filtered');
> end
> else
> %% bandpass filter for Noise cancelation of other sampling
> frequencies(Filtering)
> f1=5; %cuttoff low frequency to get rid of baseline wander
> f2=15; %cuttoff frequency to discard high frequency noise
> Wn=[f1 f2]*2/fs; % cutt off based on fs
> N = 3; % order of 3 less processing
> [a,b] = butter(N,Wn); %bandpass filtering
> ecg_h = filtfilt(a,b,ecg);
> ecg_h = ecg_h/ max( abs(ecg_h));
> if gr
> ax(3)=subplot(323);plot(ecg_h);axis tight;title('Band Pass Filtered');
> end
> end
> %% derivative filter H(z) = (1/8T)(-z^(-2) - 2z^(-1) + 2z + z^(2))
> h_d = [-1 -2 0 2 1]*(1/8);%1/8*fs
> ecg_d = conv (ecg_h ,h_d);
> ecg_d = ecg_d/max(ecg_d);
> delay = delay + 2; % delay of derivative filter 2 samples
> if gr
> ax(4)=subplot(324);plot(ecg_d);axis tight;title('Filtered with the
> derivative filter');
> end
> %% Squaring nonlinearly enhance the dominant peaks
> ecg_s = ecg_d.^2;
> if gr
> ax(5)=subplot(325);plot(ecg_s);axis tight;title('Squared');
> end
>
>
>
> %% Moving average Y(nt) = (1/N)[x(nT-(N - 1)T)+ x(nT - (N - 2)T)+...+x(nT)]
> ecg_m = conv(ecg_s ,ones(1 ,round(0.150*fs))/round(0.150*fs));
> delay = delay + 15;
>
> if gr
> ax(6)=subplot(326);plot(ecg_m);axis tight;title('Averaged with 30 samples
> length,Black noise,Green Adaptive Threshold,RED Sig Level,Red circles QRS
> adaptive threshold');
> axis tight;
> end
>
> %% Fiducial Mark
> % Note : a minimum distance of 40 samples is considered between each R wave
> % since in physiological point of view no RR wave can occur in less than
> % 200 msec distance
> [pks,locs] = findpeaks(ecg_m,'MINPEAKDISTANCE',round(0.2*fs));
>
>
>
>
> %% initialize the training phase (2 seconds of the signal) to determine the
> THR_SIG and THR_NOISE
> THR_SIG = max(ecg_m(1:2*fs))*1/3; % 0.25 of the max amplitude
> THR_NOISE = mean(ecg_m(1:2*fs))*1/2; % 0.5 of the mean signal is considered
> to be noise
> SIG_LEV= THR_SIG;
> NOISE_LEV = THR_NOISE;
>
>
> %% Initialize bandpath filter threshold(2 seconds of the bandpass signal)
> THR_SIG1 = max(ecg_h(1:2*fs))*1/3; % 0.25 of the max amplitude
> THR_NOISE1 = mean(ecg_h(1:2*fs))*1/2; %
> SIG_LEV1 = THR_SIG1; % Signal level in Bandpassed filter
> NOISE_LEV1 = THR_NOISE1; % Noise level in Bandpassed filter
> %% Thresholding and online desicion rule
>
> for i = 1 : length(pks)
>
> %% locate the corresponding peak in the filtered signal
> if locs(i)-round(0.150*fs)>= 1 && locs(i)<= length(ecg_h)
> [y_i x_i] = max(ecg_h(locs(i)-round(0.150*fs):locs(i)));
> else
> if i == 1
> [y_i x_i] = max(ecg_h(1:locs(i)));
> ser_back = 1;
> elseif locs(i)>= length(ecg_h)
> [y_i x_i] = max(ecg_h(locs(i)-round(0.150*fs):end));
> end
>
> end
>
>
> %% update the heart_rate (Two heart rate means one the moste recent and
> the other selected)
> if length(qrs_c) >= 9
>
> diffRR = diff(qrs_i(end-8:end)); %calculate RR interval
> mean_RR = mean(diffRR); % calculate the mean of 8 previous R waves
> interval
> comp =qrs_i(end)-qrs_i(end-1); %latest RR
> if comp <= 0.92*mean_RR || comp >= 1.16*mean_RR
> % lower down thresholds to detect better in MVI
> THR_SIG = 0.5*(THR_SIG);
> %THR_NOISE = 0.5*(THR_SIG);
> % lower down thresholds to detect better in Bandpass
> filtered
> THR_SIG1 = 0.5*(THR_SIG1);
> %THR_NOISE1 = 0.5*(THR_SIG1);
>
> else
> m_selected_RR = mean_RR; %the latest regular beats mean
> end
>
> end
>
> %% calculate the mean of the last 8 R waves to make sure that QRS is
> not
> % missing(If no R detected , trigger a search back) 1.66*mean
>
> if m_selected_RR
> test_m = m_selected_RR; %if the regular RR availabe use it
> elseif mean_RR && m_selected_RR == 0
> test_m = mean_RR;
> else
> test_m = 0;
> end
>
> if test_m
> if (locs(i) - qrs_i(end)) >= round(1.66*test_m)% it shows a QRS
> is missed
> [pks_temp,locs_temp] = max(ecg_m(qrs_i(end)+
> round(0.200*fs):locs(i)-round(0.200*fs))); % search back and locate the max
> in this interval
> locs_temp = qrs_i(end)+ round(0.200*fs) + locs_temp -1;
> %location
>
> if pks_temp > THR_NOISE
> qrs_c = [qrs_c pks_temp];
> qrs_i = [qrs_i locs_temp];
>
> % find the location in filtered sig
> if locs_temp <= length(ecg_h)
> [y_i_t x_i_t] =
> max(ecg_h(locs_temp-round(0.150*fs):locs_temp));
> else
> [y_i_t x_i_t] = max(ecg_h(locs_temp-round(0.150*fs):end));
> end
> % take care of bandpass signal threshold
> if y_i_t > THR_NOISE1
>
> qrs_i_raw = [qrs_i_raw locs_temp-round(0.150*fs)+
> (x_i_t - 1)];% save index of bandpass
> qrs_amp_raw =[qrs_amp_raw y_i_t]; %save amplitude of
> bandpass
> SIG_LEV1 = 0.25*y_i_t + 0.75*SIG_LEV1; %when found
> with the second thres
> end
>
> not_nois = 1;
> SIG_LEV = 0.25*pks_temp + 0.75*SIG_LEV ; %when found with
> the second threshold
> end
>
> else
> not_nois = 0;
>
> end
> end
>
>
>
>
> %% find noise and QRS peaks
> if pks(i) >= THR_SIG
>
> % if a QRS candidate occurs within 360ms of the previous
> QRS
> % ,the algorithm determines if its T wave or QRS
> if length(qrs_c) >= 3
> if (locs(i)-qrs_i(end)) <= round(0.3600*fs)
> Slope1 =
> mean(diff(ecg_m(locs(i)-round(0.075*fs):locs(i)))); %mean slope of the
> waveform at that position
> Slope2 =
> mean(diff(ecg_m(qrs_i(end)-round(0.075*fs):qrs_i(end)))); %mean slope of
> previous R wave
> if abs(Slope1) <= abs(0.5*(Slope2)) % slope
> less then 0.5 of previous R
> nois_c = [nois_c pks(i)];
> nois_i = [nois_i locs(i)];
> skip = 1; % T wave identification
> % adjust noise level in both filtered and
> % MVI
> NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1;
> NOISE_LEV = 0.125*pks(i) +
> 0.875*NOISE_LEV;
> else
> skip = 0;
> end
>
> end
> end
>
> if skip == 0 % skip is 1 when a T wave is detected
> qrs_c = [qrs_c pks(i)];
> qrs_i = [qrs_i locs(i)];
>
> % bandpass filter check threshold
> if y_i >= THR_SIG1
> if ser_back
> qrs_i_raw = [qrs_i_raw x_i]; % save index of
> bandpass
> else
> qrs_i_raw = [qrs_i_raw locs(i)-round(0.150*fs)+
> (x_i - 1)];% save index of bandpass
> end
> qrs_amp_raw =[qrs_amp_raw y_i];% save amplitude
> of bandpass
> SIG_LEV1 = 0.125*y_i + 0.875*SIG_LEV1;% adjust threshold for
> bandpass filtered sig
> end
>
> % adjust Signal level
> SIG_LEV = 0.125*pks(i) + 0.875*SIG_LEV ;
> end
>
>
> elseif THR_NOISE <= pks(i) && pks(i)<THR_SIG
>
> %adjust Noise level in filtered sig
> NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1;
> %adjust Noise level in MVI
> NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV;
>
>
>
> elseif pks(i) < THR_NOISE
> nois_c = [nois_c pks(i)];
> nois_i = [nois_i locs(i)];
>
> % noise level in filtered signal
> NOISE_LEV1 = 0.125*y_i + 0.875*NOISE_LEV1;
> %end
>
> %adjust Noise level in MVI
> NOISE_LEV = 0.125*pks(i) + 0.875*NOISE_LEV;
>
>
> end
>
>
>
>
>
> %% adjust the threshold with SNR
> if NOISE_LEV ~= 0 || SIG_LEV ~= 0
> THR_SIG = NOISE_LEV + 0.25*(abs(SIG_LEV - NOISE_LEV));
> THR_NOISE = 0.5*(THR_SIG);
> end
>
> % adjust the threshold with SNR for bandpassed signal
> if NOISE_LEV1 ~= 0 || SIG_LEV1 ~= 0
> THR_SIG1 = NOISE_LEV1 + 0.25*(abs(SIG_LEV1 - NOISE_LEV1));
> THR_NOISE1 = 0.5*(THR_SIG1);
> end
>
>
> % take a track of thresholds of smoothed signal
> SIGL_buf = [SIGL_buf SIG_LEV];
> NOISL_buf = [NOISL_buf NOISE_LEV];
> THRS_buf = [THRS_buf THR_SIG];
>
> % take a track of thresholds of filtered signal
> SIGL_buf1 = [SIGL_buf1 SIG_LEV1];
> NOISL_buf1 = [NOISL_buf1 NOISE_LEV1];
> THRS_buf1 = [THRS_buf1 THR_SIG1];
>
>
>
>
> skip = 0; %reset parameters
> not_nois = 0; %reset parameters
> ser_back = 0; %reset bandpass param
> end
>
> if gr
> hold on,scatter(qrs_i,qrs_c,'m');
> hold on,plot(locs,NOISL_buf,'--k','LineWidth',2);
> hold on,plot(locs,SIGL_buf,'--r','LineWidth',2);
> hold on,plot(locs,THRS_buf,'--g','LineWidth',2);
> if ax(:)
> linkaxes(ax,'x');
> zoom on;
> end
> end
>
>
>
>
> %% overlay on the signals
> if gr
> figure,az(1)=subplot(311);plot(ecg_h);title('QRS on Filtered Signal');axis
> tight;
> hold on,scatter(qrs_i_raw,qrs_amp_raw,'m');
> hold on,plot(locs,NOISL_buf1,'LineWidth',2,'Linestyle','--','color','k');
> hold on,plot(locs,SIGL_buf1,'LineWidth',2,'Linestyle','-.','color','r');
> hold on,plot(locs,THRS_buf1,'LineWidth',2,'Linestyle','-.','color','g');
> az(2)=subplot(312);plot(ecg_m);title('QRS on MVI signal and Noise
> level(black),Signal Level (red) and Adaptive Threshold(green)');axis tight;
> hold on,scatter(qrs_i,qrs_c,'m');
> hold on,plot(locs,NOISL_buf,'LineWidth',2,'Linestyle','--','color','k');
> hold on,plot(locs,SIGL_buf,'LineWidth',2,'Linestyle','-.','color','r');
> hold on,plot(locs,THRS_buf,'LineWidth',2,'Linestyle','-.','color','g');
> az(3)=subplot(313);plot(ecg-mean(ecg));title('Pulse train of the found QRS
> on ECG signal');axis tight;
> line(repmat(qrs_i_raw,[2 1]),repmat([min(ecg-mean(ecg))/2;
>
> max(ecg-mean(ecg))/2],size(qrs_i_raw)),'LineWidth',2.5,'LineStyle','-.','Color','r');
> linkaxes(az,'x');
> zoom on;
> end
> ##############
>
> Regards
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]
More information about the R-help
mailing list