How to define a for loop for given FINDPEAKS problem ?

[A1,L] = findpeaks(A(:,1));

A1(:,2) = A(L,2);

[A2,L1] = findpeaks(A1(:,1));

A2(:,2) = A1(L1,2);

[A3,L2] = findpeaks(A2(:,1));

A3(:,2) = A2(L2,2);

[Dyuman Joshi]

What exactly are you trying to do here?

[Parvesh]

Actually the matrix A comprises of two columns. 1st column contains acceleration (dependent variable) and the 2nd column is associated time-period (Variable).

The whole code isattached herewith.

Answer 1

[Dyuman Joshi]

Try this -

n=3;
B = [num2cell(A,1); cell(n,2)]
for k=2:n
    [B{k,1}, L] = findpeaks(B{k-1,1});
    B{k,2} = B{k-1,2}(L);
end

Here 1st row of B corresponds to A, 2nd to A1, 3rd to A2 and 4th to A3.

Use indexing to access the data.

[Parvesh]

It works and kudos to you.

[Parvesh]

Hi, there;s a small query related to the same problem. How short (with loops) the attached code can be summarised?

[Dyuman Joshi]

Please specify what you are trying to do, and what needs to be optimized?

On a cursory glance, I can say that you can read the 3 data files and do the subplots in a for loop.

type Chamoli1991VDC.m

%%%%%%%%%%%%%%%%%%%%%%% Ground Motion Selection %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Chamoli 1991 Earthquake - Acceleration (m/s/s);  %%%%%%%%%%%%%%

clc;
clear all;
close all;

%% Input:
h = 0.02;       % Time Step-Size
N = 1808;       % Number of Points
n = 4;           % Iteration needed to obtain Smooth Response Spectra
d1 = transpose(readmatrix('india.199110192123.abhat.dat','NumHeaderLines',6));□
d2 = transpose(readmatrix('india.199110192123.autta.dat','NumHeaderLines',6));□
d3 = transpose(readmatrix('india.199110192123.aghan.dat','NumHeaderLines',6));□

%% Output:□□
D1 = rmmissing(d1(:));
D2 = rmmissing(d2(:));□
D3 = rmmissing(d3(:));□
a1(:,1) = D1(1:N);□
a2(:,1) = D2(1:N);
a3(:,1) = D3(1:N);□□□□□
z1 = h*numel(a1);
z2 = h*numel(a2);
z3 = h*numel(a3);
t1 = transpose(h:h:z1);
t2 = transpose(h:h:z2);
t3 = transpose(h:h:z3);
a1(:,2) = t1(:);
a2(:,2) = t2(:);
a3(:,2) = t3(:);
A1 = [abs(a1)];
B1 = [num2cell(A1,1); cell(n,2)];
A2 = [abs(a2)];
B2 = [num2cell(A2,1); cell(n,2)];
A3 = [abs(a3)];
B3 = [num2cell(A3,1); cell(n,2)];
for k=2:n
    [B1{k,1}, L1] = findpeaks(B1{k-1,1});
    B1{k,2} = B1{k-1,2}(L1);
    [B2{k,1}, L2] = findpeaks(B2{k-1,1});
    B2{k,2} = B2{k-1,2}(L2);
    [B3{k,1}, L3] = findpeaks(B3{k-1,1});
    B3{k,2} = B3{k-1,2}(L3);
end
figure(1)
subplot(3,1,1)
plot(t1,a1(:,1),'r','LineWidth',1.25);
grid on;
xlabel('Time(in Seconds)');
ylabel('Acceleration (m/s/s)')
title('Acceleration Ground Motion - Chamoli Earthquake 1991:D1');

subplot(3,1,2)
plot(t2,a2(:,1),'b','LineWidth',1.25);
grid on;
xlabel('Time(in Seconds)');
ylabel('Acceleration (m/s/s)')
title('Acceleration Ground Motion - Chamoli Earthquake 1991:D2');

subplot(3,1,3)
plot(t3,a3(:,1),'k','LineWidth',1.25);
grid on;
xlabel('Time(in Seconds)');
ylabel('Acceleration (m/s/s)')
title('Acceleration Ground Motion - Chamoli Earthquake 1991:D3');


figure(2)
subplot(3,4,1)
plot(B1{1,2},smoothdata(B1{1,1}),'r','LineWidth',1.5);
grid on;
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')

subplot(3,4,2)
plot(B1{2,2},smoothdata(B1{2,1}),'r','LineWidth',1.5);
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')
title('Spectral Acceleration Response Spectra - Chamoli Earthquake 1991:D1')

subplot(3,4,3)
plot(B1{3,2},smoothdata(B1{3,1}),'r','LineWidth',1.5);
grid on;
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')

subplot(3,4,4)
plot(B1{4,2},smoothdata(B1{4,1}),'r','LineWidth',1.5);
grid on;
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')

subplot(3,4,5)
plot(B2{1,2},smoothdata(B2{1,1}),'b','LineWidth',1.5);
grid on;
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')

subplot(3,4,6)
plot(B2{2,2},smoothdata(B2{2,1}),'b','LineWidth',1.5);
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')
title('Spectral Acceleration Response Spectra - Chamoli Earthquake 1991:D2')

subplot(3,4,7)
plot(B2{3,2},smoothdata(B2{3,1}),'b','LineWidth',1.5);
grid on;
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')

subplot(3,4,8)
plot(B2{4,2},smoothdata(B2{4,1}),'b','LineWidth',1.5);
grid on;
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')

subplot(3,4,9)
plot(B3{1,2},smoothdata(B3{1,1}),'k','LineWidth',1.5);
grid on;
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')

subplot(3,4,10)
plot(B3{2,2},smoothdata(B3{2,1}),'k','LineWidth',1.5);
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')
title('Spectral Acceleration Response Spectra - Chamoli Earthquake 1991:D3')

subplot(3,4,11)
plot(B3{3,2},smoothdata(B3{3,1}),'k','LineWidth',1.5);
grid on;
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')

subplot(3,4,12)
plot(B3{4,2},smoothdata(B3{4,1}),'k','LineWidth',1.5);
grid on;
xlabel('Time(in Seconds)');
ylabel('Spectral Acceleration (m/s/s)')

[Parvesh]

Here, I'm reading 3 text files (these files - d1 d2 is not limited to 3 numbers, it can goo beyond that). Based on that, I'm preparing Spectral acceleration vs time curve.

Basically, these 3 files contains acceleration vs time ordinates (both + & -). A basic earthquake ground motion record and the code defined is converting them into spectral acceleration ordinates (peak absolute acceleration) and then the same are plotted down.

figure (1) comprises of comparison of 3 ground motion records while figure(2) contains comparison of spectral accelation generated from the for loop you suggested earlier

[Dyuman Joshi]

Dynamically naming variables is not a good coding practice.

Use dir to get the files in a given folder and read them in a loop -

%%%%%%%%%%% Chamoli 1991 Earthquake - Acceleration (m/s/s);  %%%%%%%%%%%%%%

%% Input:
h = 0.02;       % Time Step-Size
N = 1808;       % Number of Points
n = 4;           % Iteration needed to obtain Smooth Response Spectra

data = dir('india.199110192123.*.dat');
num = numel(d);


f1 = figure(1);
f2 = figure(2);

%text for labels
xstr = 'Time(in Seconds)';
ystr = 'Acceleration (m/s/s)';

for idx=1:num
    d = transpose(readmatrix(data(idx).name,'NumHeaderLines',6));

    %% Output:
    D = rmmissing(d(:));
    a(:,1) = D(1:N);
    t = h*transpose(1:numel(a));
    a(:,2) = t;
    A = abs(a);
    B = [num2cell(A,1); cell(n,2)];
    for k=2:n
        [B{k,1}, L] = findpeaks(B{k-1,1});
        B{k,2} = B{k-1,2}(L);
    end

    subplot(f1,3,1,idx)
    plot(t,a(:,1),'k','LineWidth',1.25);
    grid on;
    xlabel(xstr);
    ylabel(ystr);
    tstr = "Acceleration Ground Motion - Chamoli Earthquake 1991:D" + idx;
    title(tstr);

    for ij=1:4
        subplot(f2,3,4,ij+4*(idx-1))
        plot(B{ij,2},smoothdata(B{ij,1}),'r','LineWidth',1.5);
        grid on;
        xlabel(xstr);
        ylabel(ystr);
        if ij==2
            title(tstr)
        end
    end
end

[Parvesh]

it is giving error, "Unrecognized function or variable 'd'.

Error in Chamoli1991VDC02 (line 6)

num = numel(d);"

Let me share you the .dat file so that you can recheck from your end. By the way, thanks a lot man!!

[Dyuman Joshi]

Ah, that's a typo. My bad.

To correct it, replace

num = numel(d);

with

num = numel(data);