Incompatible Indices when updating matrix in the loop

Hey everyone!

I am doing some ECG analysis. I am trying to create a 2D-array that is updated for every loop iteration, but I get the following error when I run the code "Unable to perform assignment because the indices on the left side are not compatible with the size of the right side." Any idea on how to fix it?

Note: I know it is bad practice to not preallocate memory for matrices, but I am dealing with very large amounts of data and I cannot find its size.

%% Frequency Domain Analysis
fs3 = 4; % 4 hertz sampling frequency
xdft = [];
PowerSpectralDensity = [];
for i = 0: 60: (60-1)*240
    starttime2 = i;
    finishtime2 = starttime2 + 60;
    StartElement2 = starttime2 *fs+1;
    FinishElement2 = finishtime2 * fs;
    minuteSubset2 = ECG (StartElement2 : FinishElement2);
    t3 = starttime2: 1/fs3: finishtime2 - 1/fs3;

    %resampling the data to frequency = 4
    minuteSubset2_rs = resample(minuteSubset2, fs3, fs);

    N = length(minuteSubset2_rs);
    xdft = fft(minuteSubset2_rs');
    PowerSpectralDensity = (1/(fs3*N))* abs(xdft).^2;
    PowerSpectralDensity = PowerSpectralDensity(1:N/2);
    for j = 1: 240
    PDarray(i/60 +1, j) = xdft; %error in this line
    end
end

Answer 1

[William Rose]

@Mouath Abu-Daoud,

On the line in question, you are assigning vector xdft (right hand side) to one element of array PDarray(). That won't work. You must assign a vector to a row or column of an array, and the sizes must match.

Please provide values for fs and array ECG(), if you want others to be able to run this code.

[Mouath Abu-Daoud]

Thank you for your answer. How can I fix this though?

Here's my full code in case anyone needs it.

%This program analyzes ECG data over a period of 4 hours to detect if apnea
%was present

%% Opening file and extracting data
file = fopen('a01.dat');
ECG = fread(file,'int16');

fs = 100; %sampling frequency of orginal data
fs2 = 500; %new sampling frequency
t = 1/fs: 1/fs: 60;
t = t.';

% Take only the first 4 hours
ECG = ECG(1:240*60*fs)./200;
time = 1/fs:1/fs:240*60;

%% Processing and plotting data for the first minute
firstMinuteSubset = ECG (0*fs+1:60*fs);
firstMinuteSubset_rs = resample(firstMinuteSubset, fs2,fs);

    firstMinuteT1= 0: 1/fs: 60 -1/fs;
    firstMinuteT2 = 0: 1/fs2: 60 - 1/fs2;
%resampling the data to frequency = 500
    firstmMinuteSubset_rs = resample(firstMinuteSubset, fs2,fs);
%filtering the data
    d = designfilt('lowpassfir', 'filterOrder', 5, 'cutOffFrequency', 10, 'SampleRate', fs2);
    filteredFirstMinute = filter(d, firstmMinuteSubset_rs);

figure;
hold on;
plot(firstMinuteT2, filteredFirstMinute, 'r');
plot(firstMinuteT2, firstmMinuteSubset_rs, 'b');
title('Filtered and Unfiltered ECG Data');
xlabel('Time (s)');
ylabel('ECG');
%% Processing 4 hours of data
averageHRmatrix = zeros(1, 240); %preallocated matrix as in 4 hours there are 240 minutes
minimumHRmatrix = zeros(1,240);
maximumHRmatrix = zeros(1,240);
HeartRateMatrix = [];
for i = 0: 60: (60-1)*240
    starttime = i;
    finishtime = starttime + 60;
    StartElement = starttime *fs+1;
    FinishElement = finishtime * fs;
    minuteSubset = ECG (StartElement : FinishElement);
    t1= starttime: 1/fs: finishtime -1/fs;
    t2 = starttime: 1/fs2: finishtime - 1/fs2;

    %resampling the data to frequency = 500
    minuteSubset_rs = resample(minuteSubset, fs2,fs);

    %filtering the data
    d = designfilt('lowpassfir', 'filterOrder', 5, 'cutOffFrequency', 10, 'SampleRate', fs2);
    filteredECG = filter(d, minuteSubset_rs);
   
    %finding the length of RR peaks
    RRpeaks = [];
    for i2=2:length(minuteSubset) - 1 
       if (minuteSubset(i2)>1 && (minuteSubset(i2)>minuteSubset(i2+1)) && (minuteSubset(i2)>minuteSubset(i2-1)))
          RRpeaks = [RRpeaks, i2];
       end
    end
   
   
    %finding the Heart Rate

    RR_array_60s = diff(RRpeaks)/fs; 
    HR_array_60s = 1./ RR_array_60s.*60;

    HR_2 = mean(HR_array_60s)/60;

   
    HeartRateMatrix(i/60+1) = HR_2;
    minimumHRmatrix(i/60+1) = min(HR_2);
    maximumHRmatrix(i/60+1) = max(HR_2);
end
figure;
plot(HeartRateMatrix);
xlabel('Time (min)');
ylabel('Average Heart Rate (bpm)');

%% Frequency Domain Analysis
fs3 = 4; % 4 hertz sampling frequency
xdft = [];
PowerSpectralDensity = [];
for i = 0: 60: (60-1)*240
    starttime2 = i;
    finishtime2 = starttime2 + 60;
    StartElement2 = starttime2 *fs+1;
    FinishElement2 = finishtime2 * fs;
    minuteSubset2 = ECG (StartElement2 : FinishElement2);
    t3 = starttime2: 1/fs3: finishtime2 - 1/fs3;

    %resampling the data to frequency = 4
    minuteSubset2_rs = resample(minuteSubset2, fs3, fs);

    N = length(minuteSubset2_rs);
    xdft = fft(minuteSubset2_rs');
    PowerSpectralDensity = (1/(fs3*N))* abs(xdft).^2;
    PowerSpectralDensity = PowerSpectralDensity(1:N/2);
    for j = 1: 240
    PDarray(i/60 +1, j) = xdft;
    end
end

[William Rose]

@Mouath Abu-Daoud,

I think you are trying to compute the spectrum of the ECG for each one-minute interval, for four hours. I think you want array PDarray() to contain the power spectral density of the successive minutes. If that is correct, then the following should work:

PDarray(i/60+1,:)=PowerSpectralDensity;

This will assign each PSD to a new row of PDarray. I recommend pre-allocating PDarray before entering the for loop:

PDarray=zeros(240*60,M);

where M=N/2 is the expected final length of array PowerSpectralDensity.

You resample the ECG at 4 Hz. This far too slow to capture the events of the ECG.

I cannot run your code because the file a01.dat is not attached.

If this answer is helpful, then please accept it.

[Mouath Abu-Daoud]

Thank you so much, this was really helpful!

Yeah I am resampling at 4 Hz as it would be useful later on for detecting apnea. If you want I will update you with the results.

[William Rose]

@Mouath Abu-Daoud, I'm glad that my answer was helpful. Good luck with your work!