How to implement sum of sines using matricies

Hi,

Please help me to understand how to programme the matlab way... how should this example be efficiently programmed...

I am modelling a set of sin waves and would like a function that can evaluate the sum of any number of sine wave over time. The sine waves are defined by amplitude and frequency in a matrix with each row corresponding to one sine component:

M = [a1,f1;  % sin 1
     a2,f2;  % sin 2
     ...
     an,fn]; % sin n

The value of the i'th sine wave at time t is then given by:

function [Yi] = fYi(i,t)
   Yi = M(i,1) .* sin (t * M(i,2));
end

The sum of all sine wave at time t is given by:

function [Y] = fY(t)
     Y = M(1,1) * sin (t * M(1,2))+
         M(2,1) * sin (t * M(2,2))+
         M(3,1) * sin (t * M(3,2))+
         M(4,1) * sin (t * M(4,2));
end

Or if t is a scalar:

function [Y] = fY(t)
     Y = sum(M(1:4,1) * sin (t * M(1:4,2)));
end

However if the previous is called with a range of time steps e.g: t = 1:3600, then the function doesn't work due to matrix dimensions.

I basically would like the function to work with an arbitrary number of sine wave, and an arbitrary number of time points. (it should work whether called with scalar or timeseries).

Is this an instance where meshgrid could help? I cant help thinking that this should be done without looping! Please help with general advice and guidance!

Answer 1

[pfb]

You could try

 Y = a*sin(f'*t);

where a and f are the row vectors for the amplitude and frequency of your sines, and t is the row vector for time. Then f'*t is a matrix with the same rows as f' and the same columns as t, and likewise sin(f'*t).

Of course if the sizes of t and f are too large, you're going to exceed the available memory. But this is not a matlab problem. It is an hardware problem. How can you ask for "arbitrary sizes"? Your hardware is limited (although 3600*4 does not look that big).

You can very well use a loop and accumulate. There is nothing bad with it.

 Y = zeros(size(t));
 for j=1:length(a)
   Y = Y+a(j)*sin(f(j)*t);
 end

[pfb]

sorry, I did a mistake in my first reply. Now it should work.

[chris tonic]

Hi, thanks for helping! I think the first doesnt work because if there are 4 sines and 10 timesteps than f'*t yields a 4x10 matrix so multiplying this by a (also 4 row vector) fails.

Someone posted about using bsxfun to do it but have now deleted post... it seems to work. Is this method preferable to looping?:

a = 1:4;
f = 1:4;
t = 1:10;

% create 4*10 matrix of sin(f*t)
sinft = sin(bsxfun(@times, f, t));

% multiply by amplitude
Y = bsxfun(@times, a, sinft);

% sum waves
Ysum = sum(Y);

[pfb]

hi,

 Y = a*sin(f'*t);

should work if t, a and f are rows, and the last two have the same length:

 (1x4)*(4*10)=(1x10)

right?

I did a mistake in my original reply, but then I corrected it.

It seems to me that the vectorized form above does exactly the same as bsxfun, perhaps in a more intuitive way.

 sin(f'*t)

should be exactly what you call sinft. Vector-matrix multiplication does the rest.

 S=a*sinft

means that

 S(k) = a(1)*sinft(1,k)+a(2)*sinft(2,k)+...

Isn't this what you are looking for?

The memory demanding part is in sinft, which can be a very large matrix.

If you use the loop, you do not have that problem. If the sizes of t and a and f are large, I suggest the loop.

Answer 2

[Andrei Bobrov]

a = 1:4;
f = 1:4;
t = 1:10;

Ysum = sin(t(:)*f(:)')*a(:);

[pfb]

ok, you can also build a (10x4) matrix and multiply it by a (4x1) vector. Not very different from what I wrote above, though.

You get the transpose of what I get.