Summation of signals while preserving random values of the phase matrix for each individual moment in time

Question

I am summing all the harmonics of the signal in the frequency range from 0.8 to 11.1 with step 0.1 when t = const. Each individual harmonic (each signal with its own separate frequency in the range) has its own random phase. Here is the program:

clear,clc
%% Data 
x = 0:1:100; 
a = 5.76; 
omega = 0.8:0.1:11.1; 
k = omega / a; 
phase = 2*pi*rand(length(omega),1);
for i = 1:length(x)
    PHASE(i,:) = phase(:,1);
end

[X,OMEGA] = meshgrid(x,omega);
K = OMEGA / a; 

t = 0; 

y = cos(OMEGA * t - K .* X + PHASE');

y_sum = sum(y);


figure(1);
plot(x, y(1:25:end,:));
title('y(x)');
grid on;

figure(2);
plot(x, y_sum);
title('y(x)');
grid on;

A random phase is generated for each individual frequency, but relative to the x coordinates it should not change, which is why I built this loop at the beginning, which builds a matrix of random phases corresponding to these conditions (is there any other way to do this without using a loop?).

The task is as follows: for each individual moment in time t (now there are several of them), do the summation of the signals, but at the same time, so that the values of the PHASE matrix are stored for each subsequent moment in time, and not updated in a new way.

At first I tried to do this through a loop, but then I realized that every single iteration of time t it updates the value of the elements of the PHASE matrix. And as a result, the graph is not correct.

As a result, I should get graphs like this: It turned out this way only if I manually register functions y and y_sum separately for each individual moment of time, which is very inconvenient (especially if I want to change the time range).

t0 = 0; t1 = 0.1; t2 = 0.5; t3 = 1;

y0 = cos(OMEGA * t0 - K .* X + PHASE');
y1 = cos(OMEGA * t1 - K .* X + PHASE');
y2 = cos(OMEGA * t2 - K .* X + PHASE');
y3 = cos(OMEGA * t3 - K .* X + PHASE');

y_sum0 = sum(y0);
y_sum1 = sum(y1);
y_sum2 = sum(y2);
y_sum3 = sum(y3);

figure(2);
plot(x, [y_sum0',y_sum1',y_sum2',y_sum3']);
legend('t0 = 0 s', 't1 = 0.1 s', 't2 = 0.5 s', 't3 = 1 s', 'Location','eastoutside');
title('y(x)');
grid on;

How can such a task be performed in a smarter way and without multiplying functions y and y_sum for individual moments of time?

Answer 1

Hopefully the example code below answers your questions.

A few pieces of advice:

1.) Personally, I NEVER call clear or clc from a script. While using those commands is all fun and games if the script is running on your own machine, as soon as you hand the script over to a co-worker, you have no idea what those commands might erase that the co-worker didn't want erased.

2.) Do not be afraid to write your code as a function instead of a script. Doing so compartmentalizes the code so that it's easier to understand. It also lets you avoid cluttering the workspace with leftover variables (which in turn will reduce your need to call clear and clc). Lastly, if you're so inclined, you can save your functions in separate files so that you can reuse them in the future without having to copy and paste the code.

3.) Sometimes, a loop really is the best way to do things, particularly when you find yourself working with more than two dimensions of data. Don't overthink it. In this case, the example code below is likely nowhere near the complexity that would cause a noticeable performance lag, even though it's using loops.

4.) The more advanced you get with your plotting, the more helpful you will find it to retain handles to the graphics objects you create. Doing so lets you do things like alter colors, create custom legend labels, etc AFTER plotting your data. Also, if you ever want to code a MATLAB app, this is pretty much the only way to deal with plots inside the app.

I wrote the example code below in Notepad++, because I don't have MATLAB installed on the computer I'm using right now. Unfortunately, that means I'm almost certain it will have one or two dumb syntactical errors that prevent it from running. It should be pretty close, however.

% define the inputs for the problem
t = [0, 0.1, 0.5, 1];
omega = 0.8:0.1:11.1; 
x = 0:1:100;
a = 5.76;

% perform the calculation - the calcY function is declared at the bottom of the script
[y, y_sum] = calcY(t, omega, x, a);


% initialize the figure windows
yFigHandle = figure();                  % create a new figure
yAxesHandle = axes(yFigHandle);         % create a new set of cartesian axes within the new figure
title(yAxesHandle,'y(x)');
yLineObjects = gobjects(numel(y),1);    % preallocate a place holder array of graphics objects where we will eventually save the plotted lines
hold(yAxesHandle,'on');

y_sumFigHandle = figure();
y_sumAxesHandle = axes(y_sumFigHandle);
title(y_sumAxesHandle,'y_sum(x)');
y_sumLineObjects = gobjects(numel(y_sum),1);
hold(y_sumAxesHandle,'on');

% plot the data
for i = 1:numel(y)
    % By the way, I'm not sure why we're only plotting every 25th element here...there should be 114 rows in each
    % element of the yCell and that number is not divisible by 25. Once again, that's only if my back-of-the-napkin
    % math is correct, since I haven't been able to run this code.
    yLineObjects(i) = plot(yAxesHandle, x, y{i}(1:25:end,:));
    y_sumLineObjects(i) = plot(y_sumAxesHandle, x, y_sum{i});
end

% finish configuring the plots
hold(yAxesHandle,'off');
grid(yAxesHandle,'on');

hold(y_sumAxesHandle,'off');
grid(y_sumAxesHandle,'on');



function [yCell, y_sumCell] = calcY(tVector, omegaVector, xVector, a)
    % [yCell, y_sumCell] = calcY(tVector, omegaVector, xVector, a)
    %
    % Performs your calculation (which I'm still not entirely certain the purpose of, but that's not really important).
    % 
    % Returns two cell arrays, yCell and y_sumCell. They each have one column, and the same number of rows as elements
    % in tVector.
    %   Each element of yCell is a matrix of doubles with number of rows = number of elements in omegaVector and number
    %      of columns = number of elements in xVector.
    %   Each element of y_sumCell is a row vector of doubles with number of elements = number of elements in xVector.


    % ALWAYS validate your function inputs. This makes it clear to anyone else looking at this code what is expected by
    % this function, and it also prevents stupid bugs where you fed the function the wrong thing.
    arguments (Input)
        tVector (1,:) double {mustBeNonempty, mustBeReal, mustBeNonnegative}
        omegaVector (1,:) double {mustBeNonempty, mustBeReal, mustBeNonnegative}
        xVector (1,:) {mustBeNonempty, mustBeReal, mustBeNonnegative}
        a (1,1) double {mustBeReal}
    end

    % Ouput validation is optional, but I find it useful for the same reasons as above.
    arguments (Output)
        yCell (:,1) cell {mustBeNonempty}
        y_sumCell (:,1) cell {mustBeNonempty}
    end


    % determine the sizes of the output arrays
    numRows = numel(omegaVector);
    numColumns = numel(xVector);
    numCells = numel(tVector);

    % configure the inputs for the calculation
    omegaMatrix = omegaVector' .* ones(1,numColumns);       % make omegaVector a column vector and copy it to every column of omegaMatrix
    k = (omegaVector/a)';
    phase = 2*pi*rand(1,numRows);

    % preallocate the cell arrays that will be returned
    yCell = cell(numCells,1);
    y_sumCell = y;

    % loop through each element of tVector and perform the calculation.
    for i = 1:numCells
        % NOTE: If you perform element-wise multiplication between an N-element column vector and an M-element row vector,
        %   the result is an NXM matrix. Similarly, I believe that if you add a row vector to a matrix that has the same
        %   number of columns, the row vector will get added to every row.
        %   The subtraction operation is the one I'm not sure about, hence why I'm using omegaMatrix instead of just
        %   omegaVector'. Once again, I don't have MATLAB installed to test this, so I'm kind of guessing.
        yCell{i} = cos(omegaMatrix * tVector(i) - k .* xVector + phase);
        y_sumCell{i} = sum(yCell{i});
    end
end