thr3ads.net - R help - [R] Conversion of Matlab code to an R code [Mar 2015]

If this information is useful, please help other people find it:
Share via:

Abhinaba Roy

2015-Mar-23 15:10 UTC

[R] 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)<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]]

Bert Gunter

2015-Mar-23 15:23 UTC

head link

[R] Conversion of Matlab code to an R code

1. Yes. Both are Turing complete .

2. Someone here may be willing to help you if you're lucky, but the
standard recommendation is: Do your homework and learn R!  (There are
many good tutorials and resources available).

3. In future, post in plain text, not HTML, as the posting guide asks.

Cheers,
Bert

Bert Gunter
Genentech Nonclinical Biostatistics
(650) 467-7374

"Data is not information. Information is not knowledge. And knowledge
is certainly not wisdom."
Clifford Stoll




On Mon, Mar 23, 2015 at 8:10 AM, 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.

Marc Schwartz

2015-Mar-23 15:30 UTC

head link

[R] Conversion of Matlab code to an R code

> On Mar 23, 2015, at 10:10 AM, 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
<snip of 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.html

That might make your job a bit easier.

Regards,

Marc Schwartz

Roy Mendelssohn - NOAA Federal

2015-Mar-23 15:31 UTC

head link

[R] Conversion of Matlab code to an R code

On Mar 23, 2015, at 8:10 AM, 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.<snip>

If you mean something that can do an automatic conversion I am unaware of any
such program.  However, if you Google, say:

?Matlab R comparison?

you will find many links, such as ?R For Matlab Users"

http://mathesaurus.sourceforge.net/octave-r.html:

which can guide you through converting the code.

HTH,

-Roy

**********************
"The contents of this message do not reflect any position of the U.S.
Government or NOAA."
**********************
Roy Mendelssohn
Supervisory Operations Research Analyst
NOAA/NMFS
Environmental Research Division
Southwest Fisheries Science Center
***Note new address and phone***
110 Shaffer Road
Santa Cruz, CA 95060
Phone: (831)-420-3666
Fax: (831) 420-3980
e-mail: Roy.Mendelssohn at noaa.gov www: http://www.pfeg.noaa.gov/

"Old age and treachery will overcome youth and skill."
"From those who have been given much, much will be expected" 
"the arc of the moral universe is long, but it bends toward justice"
-MLK Jr.

John C Frain

2015-Mar-23 20:28 UTC

head link

[R] Conversion of Matlab code to an R code

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]]

R help - Mar 2015 - Conversion of Matlab code to an R code

[R] Conversion of Matlab code to an R code

[R] Conversion of Matlab code to an R code

[R] Conversion of Matlab code to an R code

[R] Conversion of Matlab code to an R code

[R] Conversion of Matlab code to an R code