CODEWITHSUNDEEP
matlaboctave
fearless_fool
fearless_fool

Reputation: 35189

Is there a cleaner way to generate a sum of sinusoids?

I have written a simple matlab / octave function to create the sum of sinusoids with independent amplitude, frequency and phase for each component. Is there a cleaner way to write this?

## Create a sum of cosines with independent amplitude, frequency and
## phase for each component:
##   samples(t) = SUM(A[i] * sin(2 * pi * F[i] * t + Phi[i])
## Return samples as a column vector.
##
function signal = sum_of_cosines(A = [1.0], 
                                 F = [440], 
                                 Phi = [0.0], 
                                 duration = 1.0, 
                                 sampling_rate = 44100)
  t = (0:1/sampling_rate:(duration-1/sampling_rate));
  n = length(t);
  signal = sum(repmat(A, n, 1) .* cos(2*pi*t' * F + repmat(Phi, n, 1)), 2);
endfunction

In particular, the calls to repmat() seem a bit clunky -- is there some nifty vectorization technique waiting for me to learn?

Upvotes: 2

Views: 3066

Answers (3)

Ben Voigt
Ben Voigt

Reputation: 283713

Is this the same?

signal = cos(2*pi*t' * F + repmat(Phi, n, 1)), 2) * A';

And then maybe

signal = real(exp(j*2*pi*t'*F) * (A .* exp(j*Phi))');

If you are memory constrained, this should work nicely:

e_jtheta = exp(j * 2 * pi * F / sampling_rate);
phasor = A .* exp(j*Phi);
samples = zeros(duration,1);
for k = 1:duration
    samples(k) = real((e_jtheta .^ k) * phasor');
end

Upvotes: 4

fearless_fool
fearless_fool

Reputation: 35189

Heh. When I call any of the above vectorized versions with length(A) = 10000, octave fills up VM and grinds to a halt (or at least, a slow crawl -- I haven't had the patience to wait for it to complete.

As a result, I've fallen back to the straightforward iterative version:

function signal = sum_of_cosines(A = [1.0], 
                                 F = [440], 
                                 Phi = [0.0], 
                                 duration = 1.0, 
                                 sampling_rate = 44100)
  t = (0:1/sampling_rate:(duration-1/sampling_rate));
  n = length(t);
  signal = zeros(n, 1);
  for i=1:length(A)
    samples += A(i) * cos(2*pi*t'*F(i) + Phi(i));
  endfor
endfunction

This version works plenty fast and teaches me a lesson about trying to be 'elegant'.

P.S.: This doesn't diminish my appreciation for the answers given by @BenVoigt and @chappjc -- I've learned something useful from both!

Upvotes: 0

chappjc
chappjc

Reputation: 30579

For row vectors A, F, and Phi, so you can use bsxfun to get rid of the repmat, but it is arguably uglier looking:

signal = cos(bsxfun(@plus, 2*pi*t' * F, Phi)) * A';

Upvotes: 1

Related Questions