Optimization of phase/amplitude computation loop

2 次查看（过去 30 天）
Pim Hacking
Pim Hacking 2020-9-28
评论： Pim Hacking 2020-10-7
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!

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics and Optimization 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by