Resample function is not working properly and damaging the signal

Hi,

I have a question regarding the resample function (It is the documentation under "Resample Multichannel Signal"). In this function, I realized that the function damages the signal on both ends of the given signal. Here I am attaching the code to generate 100-sample sinusoidal signal.

p = 3;
q = 2;

tx = 0:p:300-p;

x = cos(2*pi*tx./(1:5)'/100);

plot(tx,x,'.:')
title('Original')
ylim([-1.5 1.5])

However, once you resample with the following give code, the signals are simply damaged at both ends. You can see the effects on both time points close to 0 and close to 300.

ty = 0:q:300-q;

y = resample(x,p,q,'Dimension',2);

plot(ty,y,'.:')
title('Upsampled')

Therefore, could you please help me to solve this issue? Because, I am currently trying to resample my multichannel motion capture signals, which are sampled at 60 Hz or 200 Hz, in "100 Hz". However, this function will damage my signal if I use it, and this will definitely affect the outcome of my research.

I will be waiting for your help and appreciate your assistance.

Thanks!

[Walter Roberson]

Perhaps you should be using filtfilt

[goksu avdan]

So, is this "filtfilt" function eliminating the distortion of the signals for both ends?

[Walter Roberson]

I do not know; you would have to test it.

y = filtfilt(b,a,x) performs zero-phase digital filtering by processing the input data x in both the forward and reverse directions. After filtering the data in the forward direction, the function reverses the filtered sequence and runs it back through the filter. The result has these characteristics:

• Zero phase distortion.
• A filter transfer function equal to the squared magnitude of the original filter transfer function.
• A filter order that is double the order of the filter specified by b and a.

filtfilt minimizes start-up and ending transients by matching initial conditions. Do not use filtfilt with differentiator and Hilbert FIR filters, because the operation of these filters depends heavily on their phase response.

[goksu avdan]

This still affect the nature of the signal I have. It will work if I am not doing something not related about my research but I should stick with the same trend I have after resampling.

Thanks for the answer.

[Paul]

The doc page resample says that the function "compensates for the delay introduced by the filter," which is the exact language used by lowpass et. al. I'm nearly certain that the filtering in lowpass is effectively the same as filtfilt, which does extend the input signal off the edges to minimize these types of edge transients. I'm suprised that resample doesn't do the same; maybe there is something particular about the filtering in resample that precludes that (assuming it is, in fact, not done as appears to be the case).

[Bob S]

Hi,

I do have an additional comment regarding the post by @goksu avdan. There is an issue with using resample() for processing blocks of data. If I'm using resample to interpolate/decimate, I will get erroneous results at the edges. In the past, to get around this, I have used a "brute-force" approach with filter() as in the following:

while dataAvailable
    % Say x is dimensioned 1 by Nsamples and I am upsampling by a factor
    % Nup and downsampling by Ndown
    % Zero fill between samples
    xi = [x;zeros(Nup-1,Nsamples)];
    xi = xi(:).'; 
    % xi should now be dimensioned 1 by (Nup*Nsamples) with zeros inbetween the original samples
    % Now low-pass filter and decimate
    [xd, zFilt] = filter(bFilt, 1, xi, zFilt);
    xd = xd(1:Ndown:end);
end

Assuming that I use the same bFilt that would come out of resamp(), the above result will match resample() if I ignore the first floor(length(bFilt)/2/Ndown) output samples. By using the above approach, I have to deal with a transient response for the first block, after that I have continuous correct results. In order to achieve that with resample, one would have to devise some kind of overlap-save approach.

The problem with the "brute force" approach is that it is inefficient and does not take advantage of all of the zeros that a proper polyphase filtering technique, as done in resample(), would handle.

What would be ideal would be something like:

[y, zFilt] = resmaple(x, p, q, zFilt);

I do realize that depending on what the ratio of p/q looks like, it might be near impossible to get an integer number of output points without "playing tricks", but I'm ignoring that right now.

Any thoughts on this? Am I missing something? Apologies for spelling errors and any code miscues.

Thanks in advance!

Answer 1

[Star Strider]

I had not noticed the ‘end effect’ problem before (although I use resample relatively frequently, and that is clearly visible in the signals in the resample documentation). Looking at the original resampled waveforms, it occurred to me that the waveforms that were zero at both ends didn’t have that probllem. I devised a single-line linear regression for another problem recently, and decided to use it here to first ‘detrend’ the signals so that they began and ended at zero, resample them, and then ‘retrend’ them to their original orientations.

This may not be perfect, since it uses resample on the detrended signals, however it does elliminate the ‘end effect’ problem.

I am not certain what other processing you’re doing, so experiment with it to see how it works with your signals in that context —

p = 3;
q = 2;

tx = 0:p:300-p;

x = cos(2*pi*tx./(1:5)'/100);

figure
plot(tx,x,'.:')
title('Original')
ylim([-1.5 1.5])
legend(compose('%d',1:5), 'Location','best')

ty = 0:q:300-q;

% y = resample(x,p,q,'Dimension',2);

Fs = 1/(ty(2)-ty(1));

LR = @(X,Y,x,y) [X(:) ones(size(X(:)))] * ([x(:) ones(size(x(:)))] \ y(:));         % Single-Line Linear Regression
for k = 1:size(x,1)
    DeTr = LR(tx,x(k,:),tx([1 end]),x(k,[1 end])).';
    [y(k,:),ty] = resample(x(k,:)-DeTr,tx,Fs,'Dimension',2);
    ReTr = LR(ty,y(k,:),tx([end 1]),x(k,[end 1])).';
    y(k,:) = y(k,:)+ReTr;
end

figure
plot(ty,y,'.:')
title('Upsampled')
legend(compose('%d',1:5), 'Location','best')

.

[goksu avdan]

Let me try this one on my work and get back to you and update my answer.

Update: I tried to implement your model on my code; however, I got the following error. I couldn't figure it out the issue unfortunately:

"""""""

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in Extract_the_data_2 (line 53)

[y(k,:),ty] = resample(x(k,:)-DeTr,tx,Fs,'Dimension',2);

"""""""

Here is the full code; the issue is at resample function!

clear all
clc

%Change the number of subjects depending on how many subjects file you have
%ran.
n_of_subjects = 5;

for i = 3:n_of_subjects %it normally start with 1 but I only provided the sujject #3 data
    subject = ['Subject' int2str(i) '.mat']; %create the file name
    load(subject) %load the subject file (e.g., subject1, subject2...)
    [numRows, numCols] = size(s.Data); %the count of different trial numbers
    n_of_tasks = numCols;
    
    for n = 1:n_of_tasks
        result = strcmp(s.Data(n).Task, 'Walking');

        if result == 1 %trials INCLUDE walking trial
            walking = s.Data(n).Marker; %define the marker matrix
            NrNaN = sum(isnan(walking(:))) %it is able to find NaN values
            disp('true');

            %check if there is any NaN data in the trial or not!!!
            if NrNaN == 0;
                walking_trial = walking; %there is no NaN in the data, continue
                

                x = walking_trial;
                
                [numRows,numCols] = size(s.Data(n).Marker);
                
                if s.KinFreq == 200;
                    %this is for the frequency conversion
                    %from 200Hz to 100Hz
                    p = 1;
                    q = 2;
                    %number of data point multiply with p to create
                    %the timeline
                    data_count = numCols*p; 
                else
                    %this is for the frequency conversion
                    %from 60Hz to 100Hz
                    p = 5;
                    q = 3;
                    data_count = numCols*p;
                end


                tx = 0:p:data_count-p;
                ty = 0:q:data_count-q;
                
                Fs = 1/(ty(2)-ty(1));
                
                LR = @(X,Y,x,y) [X(:) ones(size(X(:)))] * ([x(:) ones(size(x(:)))] \ y(:));         % Single-Line Linear Regression
                for k = 1:size(x,1)
                    DeTr = LR(tx,x(k,:),tx([1 end]),x(k,[1 end])).';
                    [y(k,:),ty] = resample(x(k,:)-DeTr,tx,Fs,'Dimension',2);
                    ReTr = LR(ty,y(k,:),tx([end 1]),x(k,[end 1])).';
                    y(k,:) = y(k,:)+ReTr;
                end


            else
                clear walking %there is NaN in the data, clear the variable
            end

        else %trials DO NOT INCLUDE walking trial
            disp('false');
        end

    end

end

I am uploading the data into the attachment from subject 3. First two subjects data have excedded the limit of 5MB.

[Star Strider]

O.K.

[Paul]

"the waveforms that were zero at both ends didn’t have that probllem."

Are any of the waveforms zero at both ends? They all start at 1 and only one of them looks like it ends at or close to zero?

[Star Strider]

I needed to make a couple tweaks to my code to make it work in your application. They were minor, and simply involved removing the indexing from the resample output, and saving the ‘y’ vectors after the retrending step. My code itself is essentially unchanged. I also produced a plot of the result, and it appears that the ‘end effect’ is not part of those signals, and they appear to be smooth without any of the ‘end effect’.transients.

I had to change your code slightly to make it work here, and that involved commenting-out the steps for multiple file reads. I did not otherwise change it.

Try this —

% clear all
% clc

LD = load(websave('Subject3','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1328435/Subject3.mat'))

LD = struct with fields:
    s: [1×1 struct]

s = LD.s;
Data = s.Data;
Task = Data.Task;
Marker = s.Data.Marker

Marker = 87×60

0.0198 0.0484 0.0776 0.1065 0.1348 0.1620 0.1881 0.2129 0.2366 0.2591 0.2806 0.3013 0.3216 0.3417 0.3621 0.3830 0.4044 0.4265 0.4490 0.4718 0.4949 0.5183 0.5419 0.5659 0.5903 0.6152 0.6406 0.6665 0.6928 0.7196 1.5410 1.5399 1.5396 1.5402 1.5416 1.5438 1.5466 1.5498 1.5534 1.5572 1.5612 1.5653 1.5696 1.5740 1.5785 1.5828 1.5869 1.5904 1.5932 1.5949 1.5955 1.5949 1.5932 1.5903 1.5863 1.5815 1.5759 1.5699 1.5635 1.5572 0.2093 0.2219 0.2323 0.2404 0.2466 0.2515 0.2555 0.2592 0.2629 0.2667 0.2706 0.2743 0.2775 0.2800 0.2817 0.2824 0.2822 0.2814 0.2801 0.2786 0.2769 0.2752 0.2734 0.2715 0.2694 0.2671 0.2645 0.2616 0.2583 0.2548 NaN NaN NaN 0.0840 0.1095 0.1348 0.1600 0.1848 0.2094 0.2336 0.2574 0.2809 0.3040 0.3267 0.3491 0.3712 0.3931 0.4148 0.4365 0.4583 0.4804 0.5029 0.5259 0.5493 0.5733 0.5977 0.6225 0.6478 0.6735 0.6998 NaN NaN NaN 1.5419 1.5442 1.5467 1.5497 1.5531 1.5571 1.5616 1.5664 1.5715 1.5766 1.5817 1.5864 1.5907 1.5945 1.5975 1.5998 1.6010 1.6013 1.6004 1.5985 1.5954 1.5913 1.5863 1.5805 1.5742 1.5676 1.5611 NaN NaN NaN 0.0898 0.0944 0.0991 0.1036 0.1081 0.1123 0.1161 0.1194 0.1221 0.1242 0.1256 0.1263 0.1264 0.1260 0.1252 0.1241 0.1228 0.1213 0.1197 0.1179 0.1160 0.1139 0.1116 0.1090 0.1062 0.1030 0.0996 -0.2304 -0.2067 -0.1823 -0.1574 -0.1318 -0.1060 -0.0801 -0.0544 -0.0292 -0.0045 0.0195 0.0429 0.0658 0.0884 0.1107 0.1330 0.1552 0.1775 0.1998 0.2223 0.2450 0.2680 0.2913 0.3151 0.3392 0.3638 0.3888 0.4142 0.4399 0.4661 1.5015 1.4982 1.4957 1.4941 1.4934 1.4938 1.4949 1.4968 1.4994 1.5024 1.5058 1.5096 1.5136 1.5177 1.5217 1.5255 1.5288 1.5315 1.5334 1.5344 1.5344 1.5336 1.5318 1.5293 1.5260 1.5221 1.5175 1.5124 1.5071 1.5017 0.1436 0.1466 0.1504 0.1549 0.1601 0.1658 0.1717 0.1775 0.1829 0.1876 0.1915 0.1945 0.1968 0.1984 0.1995 0.2004 0.2010 0.2015 0.2018 0.2019 0.2017 0.2014 0.2008 0.2001 0.1992 0.1981 0.1968 0.1952 0.1932 0.1906 -0.3353 -0.3106 -0.2853 -0.2595 -0.2333 -0.2069 -0.1806 -0.1544 -0.1287 -0.1035 -0.0789 -0.0548 -0.0312 -0.0079 0.0152 0.0380 0.0609 0.0837 0.1065 0.1294 0.1523 0.1754 0.1986 0.2220 0.2457 0.2697 0.2941 0.3189 0.3441 0.3698

LD = load(websave('Subject3','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1328435/Subject3.mat'))

LD = struct with fields:
    s: [1×1 struct]

s = LD.s;
Data = s.Data;
Task = Data.Task;
Marker = s.Data.Marker

Marker = 87×60

0.0198 0.0484 0.0776 0.1065 0.1348 0.1620 0.1881 0.2129 0.2366 0.2591 0.2806 0.3013 0.3216 0.3417 0.3621 0.3830 0.4044 0.4265 0.4490 0.4718 0.4949 0.5183 0.5419 0.5659 0.5903 0.6152 0.6406 0.6665 0.6928 0.7196 1.5410 1.5399 1.5396 1.5402 1.5416 1.5438 1.5466 1.5498 1.5534 1.5572 1.5612 1.5653 1.5696 1.5740 1.5785 1.5828 1.5869 1.5904 1.5932 1.5949 1.5955 1.5949 1.5932 1.5903 1.5863 1.5815 1.5759 1.5699 1.5635 1.5572 0.2093 0.2219 0.2323 0.2404 0.2466 0.2515 0.2555 0.2592 0.2629 0.2667 0.2706 0.2743 0.2775 0.2800 0.2817 0.2824 0.2822 0.2814 0.2801 0.2786 0.2769 0.2752 0.2734 0.2715 0.2694 0.2671 0.2645 0.2616 0.2583 0.2548 NaN NaN NaN 0.0840 0.1095 0.1348 0.1600 0.1848 0.2094 0.2336 0.2574 0.2809 0.3040 0.3267 0.3491 0.3712 0.3931 0.4148 0.4365 0.4583 0.4804 0.5029 0.5259 0.5493 0.5733 0.5977 0.6225 0.6478 0.6735 0.6998 NaN NaN NaN 1.5419 1.5442 1.5467 1.5497 1.5531 1.5571 1.5616 1.5664 1.5715 1.5766 1.5817 1.5864 1.5907 1.5945 1.5975 1.5998 1.6010 1.6013 1.6004 1.5985 1.5954 1.5913 1.5863 1.5805 1.5742 1.5676 1.5611 NaN NaN NaN 0.0898 0.0944 0.0991 0.1036 0.1081 0.1123 0.1161 0.1194 0.1221 0.1242 0.1256 0.1263 0.1264 0.1260 0.1252 0.1241 0.1228 0.1213 0.1197 0.1179 0.1160 0.1139 0.1116 0.1090 0.1062 0.1030 0.0996 -0.2304 -0.2067 -0.1823 -0.1574 -0.1318 -0.1060 -0.0801 -0.0544 -0.0292 -0.0045 0.0195 0.0429 0.0658 0.0884 0.1107 0.1330 0.1552 0.1775 0.1998 0.2223 0.2450 0.2680 0.2913 0.3151 0.3392 0.3638 0.3888 0.4142 0.4399 0.4661 1.5015 1.4982 1.4957 1.4941 1.4934 1.4938 1.4949 1.4968 1.4994 1.5024 1.5058 1.5096 1.5136 1.5177 1.5217 1.5255 1.5288 1.5315 1.5334 1.5344 1.5344 1.5336 1.5318 1.5293 1.5260 1.5221 1.5175 1.5124 1.5071 1.5017 0.1436 0.1466 0.1504 0.1549 0.1601 0.1658 0.1717 0.1775 0.1829 0.1876 0.1915 0.1945 0.1968 0.1984 0.1995 0.2004 0.2010 0.2015 0.2018 0.2019 0.2017 0.2014 0.2008 0.2001 0.1992 0.1981 0.1968 0.1952 0.1932 0.1906 -0.3353 -0.3106 -0.2853 -0.2595 -0.2333 -0.2069 -0.1806 -0.1544 -0.1287 -0.1035 -0.0789 -0.0548 -0.0312 -0.0079 0.0152 0.0380 0.0609 0.0837 0.1065 0.1294 0.1523 0.1754 0.1986 0.2220 0.2457 0.2697 0.2941 0.3189 0.3441 0.3698

% return

%Change the number of subjects depending on how many subjects file you have
%ran.
n_of_subjects = 5;

% for i = 3:n_of_subjects %it normally start with 1 but I only provided the sujject #3 data
%     subject = ['Subject' int2str(i) '.mat']; %create the file name
%     load(subject) %load the subject file (e.g., subject1, subject2...)
    [numRows, numCols] = size(s.Data); %the count of different trial numbers
    n_of_tasks = numCols;
    % 
    for n = 1:n_of_tasks
        result = strcmp(s.Data(n).Task, 'Walking');

        if result == 1 %trials INCLUDE walking trial
            walking = s.Data(n).Marker; %define the marker matrix
            NrNaN = sum(isnan(walking(:))) %it is able to find NaN values
            disp('true');

            %check if there is any NaN data in the trial or not!!!
            if NrNaN == 0;
                walking_trial = walking; %there is no NaN in the data, continue
                
                % walking

                x = walking_trial;
                
                [numRows,numCols] = size(s.Data(n).Marker);
                
                if s.KinFreq == 200;
                    %this is for the frequency conversion
                    %from 200Hz to 100Hz
                    p = 1;
                    q = 2;
                    %number of data point multiply with p to create
                    %the timeline
                    data_count = numCols*p; 
                else
                    %this is for the frequency conversion
                    %from 60Hz to 100Hz
                    p = 5;
                    q = 3;
                    data_count = numCols*p;
                end


                tx = 0:p:data_count-p;
                ty = 0:q:data_count-q;
                
                Fs = 1/(ty(2)-ty(1));
                n
                SzX = size(x)

                LR = @(X,Y,x,y) [X(:) ones(size(X(:)))] * ([x(:) ones(size(x(:)))] \ y(:));         % Single-Line Linear Regression
                for k = 1:size(x,1)
                    DeTr = LR(tx,x(k,:),tx([1 end]),x(k,[1 end])).';
                    [y,ty] = resample(x(k,:)-DeTr,tx,Fs,'Dimension',2);
                    ReTr = LR(ty,y,tx([end 1]),x(k,[end 1])).';
                    y(k,:) = y+ReTr;
                end


            else
                clear walking %there is NaN in the data, clear the variable
            end

        else %trials DO NOT INCLUDE walking trial
            disp('false');
        end

    end

NrNaN = 9

true

NrNaN = 0

true

n = 2

SzX = 1×2

87 69

NrNaN = 0

true

n = 3

SzX = 1×2

87 76

NrNaN = 0

true

n = 4

SzX = 1×2

87 67

NrNaN = 0

true

n = 5

SzX = 1×2

87 64

NrNaN = 0

true

n = 6

SzX = 1×2

87 61

NrNaN = 0

true

n = 7

SzX = 1×2

87 77

false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false

% end

y

y = 87×129

-0.0000 -0.0003 -0.0006 -0.0010 -0.0013 -0.0016 -0.0018 -0.0020 -0.0021 -0.0021 -0.0019 -0.0016 -0.0012 -0.0006 0.0001 0.0009 0.0018 0.0026 0.0034 0.0042 0.0048 0.0052 0.0055 0.0056 0.0054 0.0051 0.0046 0.0039 0.0032 0.0024 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

figure
plot(ty,y)
grid

If you have any further problems, please do not hesitate to bring them to my attention.

.

[goksu avdan]

Hi,

I ran the code and worked without the error. However, the output matrix "y" just leaves the first and the last row of the matrix y. Rest of the data only includes Zeros. Also, the first row of the y matrix starts from the zero and ended with zero, which is the detrended version. Or am I doing something wrong?

Thanks again for your time.

[goksu avdan]

Sorry for my misexplanation. You are right, those vectors are calculated when I check the for loop step by step. However, it is deleted later on. I am currently checking the issue as well.

[Star Strider]

There is some sort of indexing glitch, however this solves it —

clear all
% clc

LD = load(websave('Subject3','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1328435/Subject3.mat'))

LD = struct with fields:
    s: [1×1 struct]

s = LD.s;
Data = s.Data;
Task = Data.Task;
Marker = s.Data.Marker

Marker = 87×60

0.0198 0.0484 0.0776 0.1065 0.1348 0.1620 0.1881 0.2129 0.2366 0.2591 0.2806 0.3013 0.3216 0.3417 0.3621 0.3830 0.4044 0.4265 0.4490 0.4718 0.4949 0.5183 0.5419 0.5659 0.5903 0.6152 0.6406 0.6665 0.6928 0.7196 1.5410 1.5399 1.5396 1.5402 1.5416 1.5438 1.5466 1.5498 1.5534 1.5572 1.5612 1.5653 1.5696 1.5740 1.5785 1.5828 1.5869 1.5904 1.5932 1.5949 1.5955 1.5949 1.5932 1.5903 1.5863 1.5815 1.5759 1.5699 1.5635 1.5572 0.2093 0.2219 0.2323 0.2404 0.2466 0.2515 0.2555 0.2592 0.2629 0.2667 0.2706 0.2743 0.2775 0.2800 0.2817 0.2824 0.2822 0.2814 0.2801 0.2786 0.2769 0.2752 0.2734 0.2715 0.2694 0.2671 0.2645 0.2616 0.2583 0.2548 NaN NaN NaN 0.0840 0.1095 0.1348 0.1600 0.1848 0.2094 0.2336 0.2574 0.2809 0.3040 0.3267 0.3491 0.3712 0.3931 0.4148 0.4365 0.4583 0.4804 0.5029 0.5259 0.5493 0.5733 0.5977 0.6225 0.6478 0.6735 0.6998 NaN NaN NaN 1.5419 1.5442 1.5467 1.5497 1.5531 1.5571 1.5616 1.5664 1.5715 1.5766 1.5817 1.5864 1.5907 1.5945 1.5975 1.5998 1.6010 1.6013 1.6004 1.5985 1.5954 1.5913 1.5863 1.5805 1.5742 1.5676 1.5611 NaN NaN NaN 0.0898 0.0944 0.0991 0.1036 0.1081 0.1123 0.1161 0.1194 0.1221 0.1242 0.1256 0.1263 0.1264 0.1260 0.1252 0.1241 0.1228 0.1213 0.1197 0.1179 0.1160 0.1139 0.1116 0.1090 0.1062 0.1030 0.0996 -0.2304 -0.2067 -0.1823 -0.1574 -0.1318 -0.1060 -0.0801 -0.0544 -0.0292 -0.0045 0.0195 0.0429 0.0658 0.0884 0.1107 0.1330 0.1552 0.1775 0.1998 0.2223 0.2450 0.2680 0.2913 0.3151 0.3392 0.3638 0.3888 0.4142 0.4399 0.4661 1.5015 1.4982 1.4957 1.4941 1.4934 1.4938 1.4949 1.4968 1.4994 1.5024 1.5058 1.5096 1.5136 1.5177 1.5217 1.5255 1.5288 1.5315 1.5334 1.5344 1.5344 1.5336 1.5318 1.5293 1.5260 1.5221 1.5175 1.5124 1.5071 1.5017 0.1436 0.1466 0.1504 0.1549 0.1601 0.1658 0.1717 0.1775 0.1829 0.1876 0.1915 0.1945 0.1968 0.1984 0.1995 0.2004 0.2010 0.2015 0.2018 0.2019 0.2017 0.2014 0.2008 0.2001 0.1992 0.1981 0.1968 0.1952 0.1932 0.1906 -0.3353 -0.3106 -0.2853 -0.2595 -0.2333 -0.2069 -0.1806 -0.1544 -0.1287 -0.1035 -0.0789 -0.0548 -0.0312 -0.0079 0.0152 0.0380 0.0609 0.0837 0.1065 0.1294 0.1523 0.1754 0.1986 0.2220 0.2457 0.2697 0.2941 0.3189 0.3441 0.3698

% return

%Change the number of subjects depending on how many subjects file you have
%ran.
n_of_subjects = 5;

% for i = 3:n_of_subjects %it normally start with 1 but I only provided the sujject #3 data
%     subject = ['Subject' int2str(i) '.mat']; %create the file name
%     load(subject) %load the subject file (e.g., subject1, subject2...)
    [numRows, numCols] = size(s.Data); %the count of different trial numbers
    n_of_tasks = numCols;
    % 
    for n = 1:n_of_tasks
        result = strcmp(s.Data(n).Task, 'Walking');

        if result == 1 %trials INCLUDE walking trial
            walking = s.Data(n).Marker; %define the marker matrix
            NrNaN = sum(isnan(walking(:))) %it is able to find NaN values
            disp('true');

            %check if there is any NaN data in the trial or not!!!
            if NrNaN == 0;
                walking_trial = walking; %there is no NaN in the data, continue
                
                % walking

                x = walking_trial;
                
                [numRows,numCols] = size(s.Data(n).Marker);
                
                if s.KinFreq == 200;
                    %this is for the frequency conversion
                    %from 200Hz to 100Hz
                    p = 1;
                    q = 2;
                    %number of data point multiply with p to create
                    %the timeline
                    data_count = numCols*p; 
                else
                    %this is for the frequency conversion
                    %from 60Hz to 100Hz
                    p = 5;
                    q = 3;
                    data_count = numCols*p;
                end


                tx = 0:p:data_count-p;
                ty = 0:q:data_count-q;
                
                Fs = 1/(ty(2)-ty(1));

                LR = @(X,Y,x,y) [X(:) ones(size(X(:)))] * ([x(:) ones(size(x(:)))] \ y(:));         % Single-Line Linear Regression
                for k = 1:size(x,1)
                    DeTr = LR(tx,x(k,:),tx([1 end]),x(k,[1 end])).';
                    [y,ty] = resample(x(k,:)-DeTr,tx,Fs,'Dimension',2);
                    ReTr = LR(ty,y,tx([end 1]),x(k,[end 1])).';
                    yc{k,:} = y+ReTr;
                    % y(k,:) = y+ReTr; 
                end


            else
                clear walking %there is NaN in the data, clear the variable
            end

        else %trials DO NOT INCLUDE walking trial
            disp('false');
        end

    end

NrNaN = 9

true

NrNaN = 0

true

NrNaN = 0

true

NrNaN = 0

true

NrNaN = 0

true

NrNaN = 0

true

NrNaN = 0

true

false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false

% end

Size_x = size(x)

Size_x = 1×2

87 77

figure
plot(tx, x)
grid
xlabel('Time')
ylabel('Signal')
title('Original Vectors')

ym = cell2mat(yc);              % Resampled Vectors
Size_ym = size(ym)

Size_ym = 1×2

87 129

figure
plot(ty,ym)
grid
xlabel('Time')
ylabel('Signal')
title('Resampled Vectors')

My code does what it is designed to do. The indexing worked correctly in my original approach, however it is getting confused here for some reason. Saving ‘y’ to a cell array and then converting that to a matrix sloves the problem.

EDIT — Re-ran code.

.

[goksu avdan]

Yes, this works. Thank you so much. I also just solved the issue. Here is the part of the code I changed, what do you think?

n
SzX = size(x)
LR = @(X,Y,x,y) [X(:) ones(size(X(:)))] * ([x(:) ones(size(x(:)))] \ y(:));         % Single-Line Linear Regression
clear a; %I clear the matrix before starting the new subject data matrix.
for k = 1:size(x,1)
    DeTr = LR(tx,x(k,:),tx([1 end]),x(k,[1 end])).';
    [y,ty] = resample(x(k,:)-DeTr,tx,Fs,'Dimension',2);

    ReTr = LR(ty,y,tx([end 1]),x(k,[end 1])).';
    a(k,:) = y+ReTr; %I assigned the new matrix name.
end

[Star Strider]

As always, my pleasure!

I completely agree! Assigning the retrended result to a different variable is definitely the key to solving it. (My problem was that I was so close to the code for so long that I failed to see that was the problem.)

I may create a function around my code and post it to the File Exchange. I need to think about that.

.

[goksu avdan]

Yes, it would be wonderful to have a function in the File Exchange.

Thanks!

[Star Strider]

Thank you!

I will work on getting it in a form I can upload. I need to add documentation as well first.

I have another one I have updated and need to test before I submit it.

.

[Stephen23]

Note that iterations of the outermost loop overwrite the data in cell array YC, leaving only the last "walking" data remaining. So for the provided data, this line:

for n = 1:n_of_tasks

is simply equivalent to

for n = 7

because all previous loops are overwritten (and the remaining loops are discarded due to NaNs). Why go to all the effort of detrending, resampling, re-trending the data for loops n = 1:6, which are all discarded?

Is that intentional?

"Yes, it would be wonderful to have a function in the File Exchange."

Using the inbuilt INTERP1 is even simpler than downloading a third-party function.

[goksu avdan]

Hi,

Thanks for the warning. I haven't done the part you mentioned. I am also going to save each of these cell arrays (or in my version it is simply assigned to the different variable) seperately, and will not be overwritten later. So, if you are asking if this is intentional, yes it was for me. Because, I didn't write that part of the code, yet. But, I have the code which can simply perform the detrending, resampling, and re-trending without any issue for each "n" in the for loop.

Regarding INTERP1, this function is not giving the same results I am looking for: This is from another post on Matlab: "The difference between interpolation (the interp1 function) and resampling (the resample function) in MATLAB is that resample is designed to resample signals, and so incorporates a FIR anti-aliasing filter. The interp1 function does not, so if you are going to do signal processing with an interpolated signal use resample, not interp1."

[Stephen23]

"Using the inbuilt INTERP1 is even simpler than downloading a third-party function."

I did not want to imply that INTERP1 replaces RESAMPLE, only that INTERP1 is a simpler replacement for Star Strider's detrend and retrend function, exactly as I showed here:

https://www.mathworks.com/matlabcentral/answers/1929455-resample-function-is-not-working-properly-and-damaging-the-signal#answer_1195665

I realize that my comment was ambiguous on this.

Answer 2

[Paul]

Hi goksu,

The observed behavior is documented at resample specifcally here: "... resample assumes that the input sequence, x, is zero before and after the samples it is given. Large deviations from zero at the endpoints of x can result in unexpected values for y." See examples at that doc page.

If you're looking for theoretically pleasing approach, resample provides means to influence transients as shown here, for example. Don't know if that approach applies to the ends of the signal or just transients internal to the signal as shown in that example.

If you're looking for a heuristic approach, one option might be to extrapolate the signal off the ends, apply resample, then chop of the extensions. I used 20 sample extensions and linear extrapolation on the example data. The actual number and method will likely depend on the actual data and the resample filter parameters. I'm only thinking about extrapolations that minimize transients at the transitions to/from the endpoints of the original signal before resampling. spline or pchip might work better. Anyway, here's what I tried:

Original signals

p = 3;
q = 2;
tx = 0:p:300-p;
x = cos(2*pi*tx./(1:5)'/100);
plot(tx,x,'.:')
title('Original')
ylim([-1.5 1.5])

Extend the signals by 20 samples via linear extrapolation.

txext = [(-20:-1)*p tx (1:20)*p+tx(end)];
xext = (interp1(tx,x.',txext,'linear','extrap')).';

Resample

y = resample(xext,p,q,'Dimension',2);
fsx = 1/p;
fsy = fsx*p/q;
ty = txext(1) + (0:(size(y,2)-1))/fsy;

Plot the extended output, note the transients have been pushed away from the central region we care about 0 < t < 300.

figure;
plot(ty,y),grid

Lop off the extensions

y(:,ty<0) = [];
ty(ty<0)  = [];
y(:,ty>tx(end)) = [];
ty(ty>tx(end))  = [];
figure
plot(ty,y),grid

[goksu avdan]

This solution looks good too. I tried another method and got some error. I am trying this one now too. Thanks.

I will update my answer once I am done with it!

[Star Strider]

If my code is not doing what you want, it is appropriate to post the signal, what my code did to it, and what the problems are with it. All I have to work with are the posted examples.

Also, I was under the impression that the ends of the signals were important data, and not to be deleted from them. That was a constraint in the approach I took.

My code detrends the signals so they are zero at the ends, resamples the detrended signals to eliminate the ‘end effect’, and then retrends the resampled signals to their original values.

[goksu avdan]

Hi Star Strider,

I posted the error under your main answer. It already is on the updated answer of mine below.

Yes, the end of data is important. You were right about it. Is Paul's answer affecting the end of the data so much or not when you extend on both ends? Because, his code worked fine for now without any error. Maybe, I can just update his answer with Spline interpolation instead of linear interpolation.

Regarding the detrended data, it is an useful approach; however, the signal I have didn't quite work well in resample function, and gave an error of "not matching dimensions". It is most likely I understand something wrong, and messed up with the code.

[Paul]

I don't think I deleted any needed data at the ends; I certainly didn't mean to if I did. The resampled, post-extraploated data is retained for all times contained within the original tx vector, including the endpoints.

Answer 3

[Stephen23]

It is easy to avoid the loop too, it is just a simple linear interpolation:

p = 3;
q = 2;
tx = (0:p:300-p).';
x = cos(2*pi*tx./(1:5)/100); % note orientation!
plot(tx,x,'.:')
title('Original')
ylim([-1.5,1.5])

ty = 0:q:300-q;
ax = interp1(tx([1,end]),x([1,end],:),tx);
ay = interp1(ty([1,end]),x([1,end],:),ty);
y = resample(x-ax,p,q)+ay;
plot(ty,y,'.:')
title('Upsampled')

[Stephen23]

Note that you can easily avoid all of those loops and all of the superfluous data processsing.

Avoid fighting MATLAB with lots of loops and IFs and the like. Think in terms of arrays and indices.

S = load('Subject3.mat');
S = S.s;
% Select valid dataset:
X = strcmp({S.Data.Task},'Walking');
Y = cellfun(@(m)any(isnan(m(:))),{S.Data.Marker});
Z = find(X&~Y,1,'last');
old = S.Data(Z).Marker;
% Sample rate:
p = interp1([60,200],[5,1],S.KinFreq);
q = interp1([60,200],[3,2],S.KinFreq);
tx = p*(0:size(old,2)-1);
Fs = 1/q;
% Resample:
DeTr = interp1(tx([1,end]),old(:,[1,end]).',tx).';
[M,ty] = resample(old-DeTr,tx,Fs,'Dimension',2);
ReTr = interp1(ty([1,end]),old(:,[1,end]).',ty).';
new = M+ReTr;
% Plot original:
figure
plot(tx, old)
grid
xlabel('Time')
ylabel('Signal')
title('Original Vectors')

% Plot resampled:
figure
plot(ty, new)
grid
xlabel('Time')
ylabel('Signal')
title('Resampled Vectors')

[goksu avdan]

Hi,

You are right, and the code you have looks a lot simplier! I can definitely avoid from the for loops and I know it. But, I did not have enough confidence and experience to perform all of these calculations. Since I do know how to do it with for loops, and I chose to do with this way. Thanks for the code!