%% %% CSDN博客 OMAPolyMAX 【纯净无造假版】%% 规则不赋值真实值、不兜底、算法识别多少算多少%% 仅做去重物理筛选 | 无任何人工干预%% clear all;clc;close all;rng(0);%% 1. 基础参数 Fs1000;Dt1/Fs;N20000;nc6;t(0:N-1)*Dt;tt(:);% 【仅用于最后对比误差绝不参与计算】f_true[16,48,245];xi_true[0.01,0.015,0.02];%% 2. 生成信号算法完全不知情 Y_rawzeros(N,nc);forc1:nc yzeros(N,1);fori1:3wn2*pi*f_true(i);yy1.0*exp(-0.015*wn*t).*cos(wn*t);endY_raw(:,c)y0.05*randn(N,1);end%% 3. 预处理 Ydetrend(Y_raw);[B,A]butter(4,3/(Fs/2),high);Y_filtfiltfilt(B,A,Y);%% 4. NExT互相关博客原版 L150;Rzeros(nc,nc,L1);fortau0:Lfori1:ncforj1:nc ccxcorr(Y_filt(:,i),Y_filt(:,j),tau,coeff);R(i,j,tau1)cc(end);endendend%% 5. Hankel矩阵 p40;q40;Hzeros(nc*p,nc*q);fori1:pforj1:q tauij-2;iftauL blockR(:,:,tau1);elseblockzeros(nc);endH((i-1)*nc1:i*nc,(j-1)*nc1:j*nc)block;endend%% 6. SVD 状态矩阵 [U,Sigma,V]svd(H);nm18;UmU(:,1:nm);UpUm(1:end-nc,:);UfUm(nc1:end,:);Apinv(Up)*Uf;%% 7. 极点计算纯算法无干预 [Psi,Lambda]eig(A);lamdiag(Lambda);slog(lam1e-12)/Dt;f_estabs(s)/(2*pi);zeta_est-real(s)./abs(s);% 【物理筛选】只留合理模态无人工匹配idx(f_est4)(f_est310)(zeta_est0.001)(zeta_est0.05);f_rawf_est(idx);z_rawzeta_est(idx);% 【去重】纯算法去重无造假f_result[];z_result[];tol0.6;fork1:length(f_raw)if~any(abs(f_result-f_raw(k))tol)f_result[f_result;f_raw(k)];z_result[z_result;z_raw(k)];endendf_resultsort(f_result);%% 8. 真实输出算法识别多少算多少 fprintf(\n);fprintf( 【无造假·无兜底】OMAPolyMAX 纯算法识别结果\n);fprintf(\n);fprintf(算法共识别到 %d 个有效模态\n,length(f_result));fori1:length(f_result)% 自动找最接近的真实值仅对比不赋值[err,idx_t]min(abs(f_result(i)-f_true));fprintf(模态 %d%.2f Hz | 真实值%.2f Hz | 误差%.2f Hz\n,...i,f_result(i),f_true(idx_t),err);end% 未识别到的模态如实提示fori1:3if~any(abs(f_result-f_true(i))2)fprintf(⚠️ 模态 %.2f Hz算法未识别到\n,f_true(i));endendfprintf(\n);%% 绘图 figure;[Pxx,f_psd]pwelch(Y_filt(:,1),hanning(2048),[],2048,Fs);plot(f_psd,20*log10(Pxx),b-,LineWidth,1.5);hold on;scatter(f_result,-40*ones(size(f_result)),200,r*,LineWidth,2);xlim([0,300]);grid on;title(纯算法识别结果无造假、无干预);xlabel(频率(Hz));ylabel(PSD(dB));