# Conversion of Matlab code to an R code Classic List Threaded 5 messages Open this post in threaded view
|

## Conversion of Matlab code to an R code

 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)
Open this post in threaded view
|

## Re: Conversion of Matlab code to an R code

Open this post in threaded view
|

## Re: Conversion of Matlab code to an R code

 In reply to this post by abhinabaroy09 > On Mar 23, 2015, at 10:10 AM, Abhinaba Roy <[hidden email]> 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 Hi, Not do the conversion automatically, certainly. I don't know that anyone will volunteer here to convert such a large volume of code, though I could be wrong of course. That being said, there are two R/Matlab references that you should leverage, if you have not already:   http://www.math.umaine.edu/~hiebeler/comp/matlabR.pdf  http://mathesaurus.sourceforge.net/octave-r.htmlThat might make your job a bit easier. Regards, Marc Schwartz ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.