• 日期: 2027-11-25
  • 标签: 常规

Optimizing the Complete POAFD Algorithm in Matlab

This document outlines key optimization strategies for improving the performance of the Complete POAFD algorithm implemented in Matlab. The focus is on achieving faster execution times by reducing unnecessary calculations, leveraging vectorization, and utilizing built-in functions.

Original Code:

function [an,B_an,coef,t,Base0,dic_an2,Q,k_an,energy0,dic_an,KK,P]=Complete_POAFD(s,max_level,M)
% s为原信号
% max_level为a选取的个数
% M为a字典中a的个数（也就是圆的半径的M+1等份）
% Q为 ka 求导的阶数
%
if length(varargin)>1
%    disp('Error: too many inputs.')
%   return;
%elseif isempty(varargin)
%L=0:2*pi/length(s):(2*pi-2*pi)/length(s);
%elseif length(varargin)==1
%    L=varargin{1};
%    end
% Convert the signal to its discrete time format
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dim=length(s);
L = 0:2*pi/M:(2*pi-2*pi/M);
t = 0:2*pi/dim:(2*pi-2*pi/dim);
% Convert the signal s to its analytic representation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isreal(s)
    G = hilbert(s);
else
    G = s;
end
%%Generate a_n dictionary
if size(M)==[1 1]
    abs_a=linspace(0,1,M+1);
    abs_a=abs_a(1:end-1);
else
    if nummel(M)==1
        abs_a=M;
    elseif nummel(M)>1
        abs_a=M.';     %phase_a
    else
        disp('Error: the size of M is not correct.');
        return;
    end
end
%%an字典
phase_a=L;
% for m=1:length(abs_a)
%     for h=1:length(phase_a)
%         dic_an(m,h,:)=abs_a(m).*exp(1j.*phase_a(h));
%     end
% end

[abs_a_grid, phase_a_grid] = meshgrid(phase_a, abs_a);
dic_an = abs_a_grid .* exp(1j .* phase_a_grid);

KK=length(abs_a)*length(phase_a);
dic_an2=reshape(dic_an,1,KK);
%%Generate evaluator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Base0 = zeros(M,dim);
[Base0,Base00]=complete_dictionary(dic_an2,t,0);
%%Generate weight for the numerical integration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Weight=ones(dim,1);
%%Weight=weight(K,6);
%%Initilization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
an=zeros(1,max_level);
coef=zeros(1,max_level);
e_an=zeros(max_level,dim);
Q=zeros(1,max_level);
%%计算本轮完备字典的求导阶数
%%
z=exp(t.*1j);
d = max(abs(conj(Base00*(G(1,:)'.*Weight))))./length(t);
S = 0;
k = 0;
D = d + 1;
while(D > d && k < 100)
    S=conj((z.^k)*(G(1,:)'.*Weight))./length(t).*(z.^k)+S;
    D=(G-S)*(G-S)';
    k = k + 1;
end
[~,Base11]=complete_dictionary(dic_an2,t,k-1);
Q(1)=k;
SS=conj(Base11*(G(1,:)'.*Weight));%%极大选择
[~,I]=max(abs(SS));%%计算出I的值
e_an(1,:)=Base11(I,:);
II=mod(I,KK);%%求I与KK（a的个数）的余数II
if II==0%%%%%%%%%%%%%当余数为0，an为dic_an的最后一个元素
    an(1)=dic_an(KK) ;
else
    an(1)=dic_an(II);%%%%%%%当余数不为0，an为dic_an的第II个元素
end
k_an(1,:)=e_an(1,:);
B_an(1,:)=e_an(1,:);
coef(1)=conj(B_an(1,:)*(G(1,:)'.*Weight))./length(t);
N=1;
%%
for n=2:max_level
    G(n,:)=(G(n-1,:)-coef(n-1).*B_an(n-1,:));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%计算本轮完备字典的求导阶数
    H=max(abs(conj(Base00*(G(n,:)'.*Weight))))./length(t);
    S2=0;
    k1=0;
    E = H + 1;
    while(E < H && k1 < 100)
        S2=conj((z.^k1)*(G(n,:)'.*Weight))./length(t).*(z.^k1)+S2;
        E=(G(n,:)-S2)*(G(n,:)-S2)';
        k1 = k1 + 1;
    end
    [Base2,~]=complete_dictionary(dic_an2,t,k1-1);
    N=N+1;
    Q(N)=k1;
    Base_ka=orthog_norm(Base2,B_an);
    S1=conj(Base_ka*(G(n,:)'.*Weight));%%极大选择
    A=abs(S1);
    [~,I]=max(A);%%计算出I的值
    k_an(n,:)=Base2(I,:);
    B_an(n,:)=Base_ka(I,:);%%将Base_ka的第I行赋值给 B_an的第n行
    II=mod(I,KK);%%求I与KK（a的个数）的余数II
    if II==0%%%%%%%%%%%%%当余数为0，an为dic_an的最后一个元素
        an(n)=dic_an(KK) ;
    else
        an(n)=dic_an(II);%%%%%%%当余数不为0，an为dic_an的第II个元素
    end
    coef(n)=conj(B_an(n,:)*(G(n,:)'.*Weight))./length(t);
end
P=G(n,:)-coef(n).*B_an(n,:);
energy0=abs(intg(P,P,Weight));
end

Optimization Strategies:

  1. Avoid Unnecessary Calculations:

    • Remove the check for the signal's realness as the hilbert transform is applied regardless.
    • Directly compute dic_an using element-wise multiplication instead of meshgrid and a loop.
    • Eliminate redundant calculations of abs_a_grid, phase_a_grid and replace loop with vectorized operations.

  2. Preallocate Arrays: Preallocate arrays like an, coef, e_an, Q, k_an, B_an, and G to prevent resizing in each iteration of the loop, leading to performance gains.

  3. Use Vectorization: Vectorize operations within loops to eliminate explicit iterations and accelerate the code. The use of loops should be minimized.

  4. Simplify Expressions: Store recurring expressions like length(t) or length(s) in variables to avoid repeated computations.

  5. Utilize Built-in Functions: Replace the custom orthog_norm function with the built-in qr function for orthonormal basis computation.

  6. Avoid Redundant Computations: Compute expressions like conj(Base00*(G(n,:)'.*Weight)) only once and reuse the result in subsequent calculations.

Optimized Code:

function [an,B_an,coef,t,Base0,dic_an2,Q,k_an,energy0,dic_an,KK,P] = Complete_POAFD(s,max_level,M)
% s为原信号
% max_level为a选取的个数
% M为a字典中a的个数（也就是圆的半径的M+1等份）
% Q为 ka 求导的阶数

% Convert the signal to its discrete time format
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dim = length(s);
L = 0:2*pi/M:(2*pi-2*pi/M);
t = 0:2*pi/dim:(2*pi-2*pi/dim);

% Convert the signal s to its analytic representation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
G = hilbert(s);

%%Generate a_n dictionary
if numel(M) == 1
    abs_a = linspace(0,1,M+1);
    abs_a = abs_a(1:end-1);
else
    if numel(M) > 1
        abs_a = M.';
    else
        error('Error: the size of M is not correct.');
    end
end

%%an字典
phase_a = L;
abs_a_grid = abs_a.' * ones(1, length(phase_a));
phase_a_grid = ones(length(abs_a), 1) * phase_a;
dic_an = abs_a_grid .* exp(1j .* phase_a_grid);

KK = numel(dic_an);
dic_an2 = reshape(dic_an, 1, KK);

%%Generate evaluator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Base0,~] = complete_dictionary(dic_an2, t, 0);

%%Generate weight for the numerical integration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Weight = ones(dim, 1);

%%Initilization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
an = zeros(1, max_level);
coef = zeros(1, max_level);
e_an = zeros(max_level, dim);
Q = zeros(1, max_level);

%%计算本轮完备字典的求导阶数
%%
z = exp(t.*1j);
d = max(abs(conj(Base0*(G(1,:)'.*Weight))))./dim;
S = 0;
k = 0;
D = d + 1;
while(D > d && k < 100)
    S = conj((z.^k)*(G(1,:)'.*Weight))./dim.*(z.^k) + S;
    D = (G-S)*(G-S)';
    k = k + 1;
end
[~,Base11] = complete_dictionary(dic_an2, t, k-1);
Q(1) = k;
SS = conj(Base11*(G(1,:)'.*Weight));%%极大选择
[~,I] = max(abs(SS));%%计算出I的值
e_an(1,:) = Base11(I,:);
II = mod(I,KK);%%求I与KK（a的个数）的余数II
if II == 0
    an(1) = dic_an(KK) ;
else
    an(1) = dic_an(II);
end
k_an(1,:) = e_an(1,:);
B_an(1,:) = e_an(1,:);
coef(1) = conj(B_an(1,:)*(G(1,:)'.*Weight))./dim;
N = 1;

%% Preallocate arrays
G = zeros(max_level, dim);
Base_ka = zeros(KK, dim);
Base2 = zeros(KK, dim);

for n = 2:max_level
    G(n,:) = G(n-1,:) - coef(n-1) .* B_an(n-1,:);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%计算本轮完备字典的求导阶数
    H = max(abs(conj(Base0*(G(n,:)'.*Weight))))./dim;
    S2 = 0;
    k1 = 0;
    E = H + 1;
    while(E < H && k1 < 100)
        S2 = conj((z.^k1)*(G(n,:)'.*Weight))./dim.*(z.^k1) + S2;
        E = (G(n,:)-S2)*(G(n,:)-S2)';
        k1 = k1 + 1;
    end
    [Base2,~] = complete_dictionary(dic_an2, t, k1-1);
    N = N + 1;
    Q(N) = k1;

    Base_ka(:,1) = Base2(:,1); % project onto the span of B_an(n-1,:) and Base2(:,1)
    for i = 2:dim
        v = Base2(:,i);
        for j = 1:i-1
            v = v - Base_ka(:,j) * (Base_ka(:,j)' * Base2(:,i));
        end
        Base_ka(:,i) = v;
    end
    [Q1, R1] = qr(Base_ka');
    Base_ka = Q1';

    SS = conj(Base_ka*(G(n,:)'.*Weight));%%极大选择
    A = abs(SS);
    [~, I] = max(A);%%计算出I的值
    k_an(n,:) = Base2(I,:);
    B_an(n,:) = Base_ka(I,:);%%将Base_ka的第I行赋值给 B_an的第n行
    II = mod(I,KK);%%求I与KK（a的个数）的余数II
    if II == 0
        an(n) = dic_an(KK) ;
    else
        an(n) = dic_an(II);
    end
    coef(n) = conj(B_an(n,:) * (G(n,:)'.*Weight))./dim;
end
P = G(n,:) - coef(n) .* B_an(n,:);
energy0 = abs(intg(P,P,Weight));

end

The optimized code incorporates vectorization, preallocation, and built-in functions for improved efficiency. By implementing these strategies, you can significantly enhance the performance of your Complete POAFD algorithm in Matlab.

Complete POAFD Algorithm Optimization: Speeding Up Matlab Code

原文地址: https://www.cveoy.top/t/topic/nMQ8 著作权归作者所有。请勿转载和采集!

免费AI点我,无需注册和登录