Optimization of phase/amplitude computation loop

Hi all,

I am stuck trying to optimize my for-loop using matrix arithmetic. Perhaps you guys can help me out with optimizing this loop and/or a different approach to compute phase/amplitude. Below you'll find my current loop (both the unoptimized / my attempt at optimization). Note that I'm trying to process a matrix of signals! These signals are stored in a KxLxN matrix (where N are the number of samples in each signal). Note that typically L >> N > K.

I am working on tooling which includes determining parameters of damped sinusoids. I already have methods to derive all parameters, I'm only unhappy with the performance of my amplitude estimation. I merely compute the phase in this way because I already spent the computational effort to compute the amplitude (and the phase is only a minor addition to this). I'm open to other suggestions, both speed and accuracy are important! At this point in the algorithm I've already computed the frequency and damping. I also have the original (noisy) and FFT signal available (note I've applied some spectral leakage windowing in the FFT!).

Original Loop

C = 2*pi*frequency / sample_freq;
B = -damping .* C;  

X = (0:(N-1))';
D = zeros(K, L, 2);
phase = zeros(K, L);
for k = 1:K
    for l = 1:L
        exp_B = exp(B(k,l)*X);
        C_X   = C(k,l)*X;
        U     = [exp_B .* sin(C_X),...
                 exp_B .* cos(C_X)]; % compute basis vectors

        D(k,l,:) = U \ reshape(original(k,l,:),[N 1]); % Solve U with Original signal
        phase(k,l) = atan2d(D(k,l,2),D(k,l,1)); % compute phase
    end
end

amplitude = sqrt(D(:,:,1).^2 + D(:,:,2).^2); % = norm(d);

So I tried to do the obvious thing, which is to get rid of the loop (as much as possible!). I got stuck with removing "solving the system" (pinv?) and atan2d from the loop, as I am not using those frequently. Things actually got slower in my optimization attempt due to the addtion of a permute().

(not so) Optimized Loop

C = 2*pi*frequency / sample_freq;
B = -damping .* C;  

X = repelem( (0:(N-1))', 1, K, L) ); 
X = permute(X, [2 3 1]);
D = zeros(K, L, 2);
phase = zeros(K, L);

exp_B = exp(B.*X);
C_X = C .* X;

U1 = exp_B .* sin(C_X);
U2 = exp_B .* cos(C_X);

for k = 1:K
    for l = 1:L
        D(k,l,:) = permute([U1(k,l,:); U2(k,l,:)],[3 1 2]) \ reshape(original(k,l,:),[N 1]);
        phase(k,l) = atan2d(D(k,l,2),D(k,l,1)); % compute phase
    end
end

amplitude = sqrt(D(:,:,1).^2 + D(:,:,2).^2); % = norm(d);

Any help is much appreciated!

[Pim Hacking]

bump