E - Dominant Indices

本文详细解析了一种基于树形结构的算法问题,旨在计算每个节点的深度数组及其主导索引。通过递归深度优先搜索(DFS)策略,结合map数据结构来记录每个节点子树中各个深度的节点数,最终确定每个节点的主导深度。代码实现采用C++,展示了如何高效处理大规模树形数据。

You are given a rooted undirected tree consisting of n vertices. Vertex 1 is the root.

Let’s denote a depth array of vertex x as an infinite sequence [dx,0,dx,1,dx,2,…], where dx,i is the number of vertices y such that both conditions hold:

x is an ancestor of y;
the simple path from x to y traverses exactly i edges.
The dominant index of a depth array of vertex x (or, shortly, the dominant index of vertex x) is an index j such that:

for every k<j, dx,k<dx,j;
for every k>j, dx,k≤dx,j.
For every vertex in the tree calculate its dominant index.

Input
The first line contains one integer n (1≤n≤106) — the number of vertices in a tree.

Then n−1 lines follow, each containing two integers x and y (1≤x,y≤n, x≠y). This line denotes an edge of the tree.

It is guaranteed that these edges form a tree.

Output
Output n numbers. i-th number should be equal to the dominant index of vertex i.

Examples
Input
4
1 2
2 3
3 4
Output
0
0
0
0
Input
4
1 2
1 3
1 4
Output
1
0
0
0
Input
4
1 2
2 3
2 4
Output
2
1
0
0
题意变化一下就是找到每个结点子树中同一深度最多的那个深度,相同时取最近的,用map记录每个结点子树结点在每个深度下的数量,用两个数组分别记录深度最多的结点数量和最近深度用来更新。

#include<bits/stdc++.h>
using namespace std;
const int N=1e6+5;
int rt[N],dep[N];
int ma[N],num[N];
int ans[N];
vector<int> edge[N];
map<int,int> m[N];
map<int,int>::iterator it;

int merge(int x,int y)
{
	if(m[x].size()<m[y].size())
	swap(x,y);
	
	for(it=m[y].begin();it!=m[y].end();it++)
	{
		m[x][it->first]+=it->second;
		if(m[x][it->first]>num[x]||(m[x][it->first]==num[x]&&(it->first)<=ma[x]))
		{
			num[x]=m[x][it->first];
			ma[x]=it->first;
		}
	}
	return x;
}

void dfs(int x,int fa)
{
	dep[x]=dep[fa]+1;
	ma[rt[x]]=dep[x];
	num[rt[x]]=1;
	m[rt[x]][dep[x]]++;
	for(int i=0;i<edge[x].size();i++)
	if(edge[x][i]!=fa)
	dfs(edge[x][i],x);
	
	for(int i=0;i<edge[x].size();i++)
	{
		if(edge[x][i]!=fa)
		rt[x]=merge(rt[x],rt[edge[x][i]]);
	}
	ans[x]=ma[rt[x]];
}

int n,a,b;
int main()
{
	scanf("%d",&n);
	dep[0]=0;
	for(int i=1;i<n;i++)
	{
		scanf("%d%d",&a,&b);
		edge[a].push_back(b);
		edge[b].push_back(a);
		rt[i]=i;
	}
	rt[n]=n;
	dfs(1,0);
	for(int i=1;i<=n;i++)
	printf("%d\n",ans[i]-dep[i]);
}
import os import re import pandas as pd import numpy as np import tensorflow as tf from sklearn.preprocessing import StandardScaler, OneHotEncoder from sklearn.compose import ColumnTransformer from sklearn.model_selection import train_test_split from tensorflow.keras.layers import Input, Dense, Concatenate, Multiply, GlobalAveragePooling1D, Conv1D from tensorflow.keras.models import Model from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau import matplotlib.pyplot as plt import seaborn as sns # --------------------------- # 1. 数据加载与完整性验证 # --------------------------- def find_gait_data_files(data_dir='gait_data'): """自动发现所有步态数据文件""" file_map = {} patterns = [ r'koa_?(\w+)\.npy', # koa_t25.npy 或 koat25.npy r'KOA_?(\w+)\.npy', # KOA_T25.npy r'patient_?(\w+)\.npy', # patient_t25.npy r'(\w+)_gait\.npy' # t25_gait.npy ] for filename in os.listdir(data_dir): if filename.endswith('.npy'): for pattern in patterns: match = re.match(pattern, filename, re.IGNORECASE) if match: patient_id = match.group(1) file_map[patient_id.lower()] = os.path.join(data_dir, filename) break return file_map def load_and_validate_data(metadata_path='patient_metadata.csv', data_dir='gait_data'): """加载并验证数据完整性""" # 加载元数据 metadata = pd.read_csv(metadata_path, encoding='utf-8-sig') # 清理患者ID metadata['clean_id'] = metadata['patient_id'].apply( lambda x: re.sub(r'[^\w]', '', str(x)).str.lower() # 获取步态文件映射 gait_file_map = find_gait_data_files(data_dir) # 添加文件路径到元数据 metadata['gait_path'] = metadata['clean_id'].map(gait_file_map) # 标记缺失数据 metadata['data_complete'] = metadata['gait_path'].notnull() # 分离完整和不完整数据 complete_data = metadata[metadata['data_complete']].copy() incomplete_data = metadata[~metadata['data_complete']].copy() print(f"总患者数: {len(metadata)}") print(f"完整数据患者: {len(complete_data)}") print(f"缺失步态数据患者: {len(incomplete_data)}") if not incomplete_data.empty: print("\n缺失步态数据的患者:") print(incomplete_data[['patient_id', 'koa_grade']]) incomplete_data.to_csv('incomplete_records.csv', index=False) return complete_data # --------------------------- # 2. 特征工程 # --------------------------- def preprocess_static_features(metadata): """预处理静态元数据特征""" # 定义特征类型 numeric_features = ['age', 'height_cm', 'weight_kg', 'bmi'] categorical_features = ['gender', 'dominant_leg', 'affected_side'] # 确保所有列存在 for col in numeric_features + categorical_features: if col not in metadata.columns: metadata[col] = np.nan # 填充缺失值 metadata[numeric_features] = metadata[numeric_features].fillna(metadata[numeric_features].median()) metadata[categorical_features] = metadata[categorical_features].fillna('unknown') # 创建预处理管道 preprocessor = ColumnTransformer( transformers=[ ('num', StandardScaler(), numeric_features), ('cat', OneHotEncoder(handle_unknown='ignore'), categorical_features) ], remainder='drop' # 只处理指定特征 ) # 应用预处理 static_features = preprocessor.fit_transform(metadata) feature_names = ( numeric_features + list(preprocessor.named_transformers_['cat'].get_feature_names_out(categorical_features)) return static_features, feature_names def extract_gait_features(gait_path): """从步态数据中提取特征""" gait_data = np.load(gait_path) # 1. 时域特征 mean_values = np.mean(gait_data, axis=0) std_values = np.std(gait_data, axis=0) max_values = np.max(gait_data, axis=0) min_values = np.min(gait_data, axis=0) # 2. 频域特征 (FFT) fft_values = np.abs(np.fft.rfft(gait_data, axis=0)) dominant_freq = np.argmax(fft_values, axis=0) spectral_energy = np.sum(fft_values**2, axis=0) # 3. 生物力学特征 # 左右腿不对称性 asymmetry_index = np.mean(np.abs(gait_data[:, :5] - gait_data[:, 5:10]), axis=0) # 关节协调性 (髋-膝相位差) hip_knee_phase_diff = np.arctan2(gait_data[:, 2], gait_data[:, 1]) - np.arctan2(gait_data[:, 7], gait_data[:, 6]) phase_diff_mean = np.mean(hip_knee_phase_diff) phase_diff_std = np.std(hip_knee_phase_diff) # 步态周期对称性 gait_cycle_symmetry = np.corrcoef(gait_data[:, :5].flatten(), gait_data[:, 5:10].flatten())[0, 1] # 组合所有特征 gait_features = np.concatenate([ mean_values, std_values, max_values, min_values, dominant_freq, spectral_energy, asymmetry_index, [phase_diff_mean, phase_diff_std, gait_cycle_symmetry] ]) # 特征名称 base_names = [ 'R.Hip_AP', 'R.Knee_Moment', 'R.KneeFlex', 'R.Gluteus', 'R.Vastus', 'L.Hip_AP', 'L.Knee_Moment', 'L.KneeFlex', 'L.Gluteus', 'L.Vastus' ] feature_names = [] for prefix in ['mean_', 'std_', 'max_', 'min_']: feature_names.extend([prefix + name for name in base_names]) feature_names.extend(['dom_freq_' + name for name in base_names]) feature_names.extend(['spec_energy_' + name for name in base_names]) feature_names.extend(['asym_' + name for name in base_names]) feature_names.extend(['phase_diff_mean', 'phase_diff_std', 'gait_symmetry']) return gait_features, feature_names def biomechanical_feature_enhancement(static_df, gait_features): """基于临床知识的特征增强""" enhanced_features = gait_features.copy() # BMI影响肌肉力量 bmi = static_df['bmi'].values[0] bmi_factor = 1 + (bmi - 25) * 0.01 # 每增加1 BMI单位增加1%肌肉力量 # 肌肉力量通道索引 (R.Gluteus, R.Vastus, L.Gluteus, L.Vastus) muscle_indices = [3, 4, 8, 9] for idx in muscle_indices: # 应用BMI调整 enhanced_features[idx] *= bmi_factor # mean enhanced_features[idx + 10] *= bmi_factor # std enhanced_features[idx + 20] *= bmi_factor # max enhanced_features[idx + 30] *= bmi_factor # min enhanced_features[idx + 50] *= bmi_factor # spectral energy # 年龄影响步态稳定性 age = static_df['age'].values[0] stability_factor = max(0.8, 1 - (age - 60) * 0.005) # 60岁以上每岁减少0.5%稳定性 # 稳定性相关特征 stability_indices = [ 1, 6, # 膝力矩 2, 7, # 膝屈曲 40, 45, 50, 55, # 不对称性和相位差特征 len(gait_features) - 3, len(gait_features) - 2, len(gait_features) - 1 # 相位和对称性 ] for idx in stability_indices: enhanced_features[idx] *= stability_factor return enhanced_features # --------------------------- # 3. 特征融合模型 # --------------------------- def create_fusion_model(static_dim, gait_dim, num_classes): """创建特征融合模型""" # 输入层 static_input = Input(shape=(static_dim,), name='static_input') gait_input = Input(shape=(gait_dim,), name='gait_input') # 静态特征处理 static_stream = Dense(64, activation='relu')(static_input) static_stream = Dense(32, activation='relu')(static_stream) # 步态特征处理 gait_stream = Dense(128, activation='relu')(gait_input) gait_stream = Dense(64, activation='relu')(gait_stream) # 注意力融合机制 attention_scores = Dense(64, activation='sigmoid')(gait_stream) attended_static = Multiply()([static_stream, attention_scores]) # 特征融合 fused = Concatenate()([attended_static, gait_stream]) # 分类层 x = Dense(128, activation='relu')(fused) x = Dense(64, activation='relu')(x) output = Dense(num_classes, activation='softmax', name='koa_grade')(x) return Model(inputs=[static_input, gait_input], outputs=output) def create_1d_cnn_fusion_model(static_dim, gait_timesteps, gait_channels, num_classes): """创建包含原始时间序列的融合模型""" # 输入层 static_input = Input(shape=(static_dim,), name='static_input') gait_input = Input(shape=(gait_timesteps, gait_channels), name='gait_input') # 静态特征处理 static_stream = Dense(64, activation='relu')(static_input) static_stream = Dense(32, activation='relu')(static_stream) # 时间序列特征提取 gait_stream = Conv1D(64, 7, activation='relu', padding='same')(gait_input) gait_stream = Conv1D(128, 5, activation='relu', padding='same')(gait_stream) gait_stream = GlobalAveragePooling1D()(gait_stream) # 特征融合 fused = Concatenate()([static_stream, gait_stream]) # 分类层 x = Dense(128, activation='relu')(fused) x = Dense(64, activation='relu')(x) output = Dense(num_classes, activation='softmax', name='koa_grade')(x) return Model(inputs=[static_input, gait_input], outputs=output) # --------------------------- # 4. 训练与评估 # --------------------------- def train_model(metadata, model_type='features'): """训练特征融合模型""" # 预处理静态特征 static_features, static_feature_names = preprocess_static_features(metadata) # 提取步态特征 gait_features_list = [] gait_feature_names = None gait_data_list = [] # 原始时间序列数据 for idx, row in metadata.iterrows(): gait_path = row['gait_path'] gait_features, gait_feature_names = extract_gait_features(gait_path) # 应用生物力学增强 enhanced_features = biomechanical_feature_enhancement(row.to_frame().T, gait_features) gait_features_list.append(enhanced_features) # 保存原始时间序列数据(用于CNN模型) gait_data = np.load(gait_path) gait_data_list.append(gait_data) gait_features = np.array(gait_features_list) # 目标变量 y = metadata['koa_grade'].values num_classes = len(np.unique(y)) # 划分训练测试集 X_static_train, X_static_test, X_gait_train, X_gait_test, y_train, y_test = train_test_split( static_features, gait_features, y, test_size=0.2, random_state=42, stratify=y ) # 创建模型 if model_type == 'features': model = create_fusion_model( static_dim=static_features.shape[1], gait_dim=gait_features.shape[1], num_classes=num_classes ) train_input = [X_static_train, X_gait_train] test_input = [X_static_test, X_gait_test] else: # 使用原始时间序列 # 填充时间序列到相同长度 max_timesteps = max(len(data) for data in gait_data_list) X_gait_padded = np.array([np.pad(data, ((0, max_timesteps - len(data)), (0, 0)), 'constant') for data in gait_data_list]) # 重新划分数据集 X_static_train, X_static_test, X_gait_train, X_gait_test, y_train, y_test = train_test_split( static_features, X_gait_padded, y, test_size=0.2, random_state=42, stratify=y ) model = create_1d_cnn_fusion_model( static_dim=static_features.shape[1], gait_timesteps=X_gait_padded.shape[1], gait_channels=X_gait_padded.shape[2], num_classes=num_classes ) train_input = [X_static_train, X_gait_train] test_input = [X_static_test, X_gait_test] # 编译模型 model.compile( optimizer='adam', loss='sparse_categorical_crossentropy', metrics=['accuracy'] ) # 训练模型 history = model.fit( train_input, y_train, validation_split=0.1, epochs=100, batch_size=16, callbacks=[ EarlyStopping(patience=20, restore_best_weights=True), ReduceLROnPlateau(factor=0.5, patience=10) ] ) # 评估模型 test_loss, test_acc = model.evaluate(test_input, y_test) print(f"\n测试准确率: {test_acc:.4f}") # 特征重要性分析 if model_type == 'features': analyze_feature_importance(model, X_static_test, X_gait_test, y_test, static_feature_names, gait_feature_names) return model, history def analyze_feature_importance(model, X_static, X_gait, y_true, static_names, gait_names): """分析特征重要性""" # 获取模型预测 y_pred = model.predict([X_static, X_gait]) y_pred_classes = np.argmax(y_pred, axis=1) # 计算混淆矩阵 cm = confusion_matrix(y_true, y_pred_classes) plt.figure(figsize=(10, 8)) sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', xticklabels=["G0", "G1", "G2", "G3", "G4"], yticklabels=["G0", "G1", "G2", "G3", "G4"]) plt.xlabel('预测分级') plt.ylabel('真实分级') plt.title('KL分级混淆矩阵') plt.savefig('confusion_matrix.png', dpi=300) plt.close() # SHAP特征重要性分析 try: import shap # 创建解释器 explainer = shap.KernelExplainer(model.predict, [X_static[:50], X_gait[:50]]) # 计算SHAP值 shap_values = explainer.shap_values([X_static[:50], X_gait[:50]]) # 组合特征名称 all_feature_names = static_names + gait_names # 可视化 plt.figure(figsize=(12, 8)) shap.summary_plot(shap_values, feature_names=all_feature_names, plot_type='bar') plt.tight_layout() plt.savefig('feature_importance.png', dpi=300) plt.close() except ImportError: print("SHAP库未安装,跳过特征重要性分析") # --------------------------- # 5. 主工作流程 # --------------------------- def main(): print("="*50) print("KOA分级诊疗平台 - 特征融合系统") print("="*50) # 步骤1: 加载并验证数据 print("\n[步骤1] 加载数据并验证完整性...") complete_data = load_and_validate_data() if len(complete_data) < 20: print("\n错误: 完整数据样本不足,无法训练模型") return # 步骤2: 训练特征融合模型 print("\n[步骤2] 训练特征融合模型...") print("选项1: 使用提取的特征 (更快)") print("选项2: 使用原始时间序列 (更准确)") choice = input("请选择模型类型 (1/2): ").strip() model_type = 'features' if choice == '1' else 'time_series' model, history = train_model(complete_data, model_type) # 保存模型 model.save(f'koa_fusion_model_{model_type}.h5') print(f"模型已保存为 koa_fusion_model_{model_type}.h5") # 绘制训练历史 plt.figure(figsize=(12, 6)) plt.subplot(1, 2, 1) plt.plot(history.history['accuracy'], label='训练准确率') plt.plot(history.history['val_accuracy'], label='验证准确率') plt.title('模型准确率') plt.ylabel('准确率') plt.xlabel('Epoch') plt.legend() plt.subplot(1, 2, 2) plt.plot(history.history['loss'], label='训练损失') plt.plot(history.history['val_loss'], label='验证损失') plt.title('模型损失') plt.ylabel('损失') plt.xlabel('Epoch') plt.legend() plt.tight_layout() plt.savefig('training_history.png', dpi=300) print("训练历史图已保存为 training_history.png") print("\n特征融合模型训练完成!") if __name__ == "__main__": main()我那个步态数据文件是Excel的
07-24
提取所有数学关系 import numpy as np import matplotlib.pyplot as plt import pandas as pd import os import pywt from scipy.signal import savgol_filter from scipy.interpolate import interp1d import sys import io import json plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False def create_result_folder(folder_name="results"): if not os.path.exists(folder_name): os.makedirs(folder_name) return folder_name def save_plot(save_path, filename, fig=None): if fig is not None: fig.savefig(os.path.join(save_path, filename), dpi=300, bbox_inches='tight') plt.close(fig) else: plt.savefig(os.path.join(save_path, filename), dpi=300, bbox_inches='tight') plt.close() path = os.path.dirname(os.path.abspath(__file__)) # 1. 数据读取与异常值剔除 def get_data(file_path): file_name = os.path.basename(file_path) data = pd.read_excel(file_path, header=0) if data.isna().any().any(): print(f"{file_name}存在空值:") print(data.isna().sum()) else: print(f"{file_name}无空值") wavenumber = data.iloc[:, 0].values reflectivity = data.iloc[:, 1].values negative_num = np.sum(reflectivity < 0) over_one_num = np.sum(reflectivity > 100) print(f"{file_name}反射率小于0的个数:{negative_num}") print(f"{file_name}反射率大于100的个数:{over_one_num}") return wavenumber, reflectivity def outlier_removal(wavenumber, reflectivity): valid_value = (reflectivity >= 0) & (reflectivity <= 100) return wavenumber[valid_value], reflectivity[valid_value] # 2. 数据平滑与去噪 def apply_savgol_filter(reflectivity, window_length=11, polyorder=3): return savgol_filter(reflectivity, window_length=window_length, polyorder=polyorder) def wavelet_denoise(signal, wavelet='db4', level=1, mode_ext='constant'): coeffs = pywt.wavedec(signal, wavelet, level=level, mode=mode_ext) threshold = np.std(coeffs[-level]) * np.sqrt(2 * np.log(len(signal))) coeffs = [pywt.threshold(c, threshold, mode='soft') for c in coeffs] reconstructed = pywt.waverec(coeffs, wavelet, mode=mode_ext) return reconstructed def equal_length(signal, target_length): x = np.linspace(0, 1, len(signal)) x_target = np.linspace(0, 1, target_length) f = interp1d(x, signal, kind='linear') return f(x_target) def normalize(reflectivity): return reflectivity / 100.0 # 3. 导数法检测波峰波谷 def find_peaks_by_derivative(wavenumber, reflectivity, threshold=1e-3): first_derivative = np.gradient(reflectivity, wavenumber) second_derivative = np.gradient(first_derivative, wavenumber) peak_indices = np.where((np.abs(first_derivative) < threshold) & (second_derivative < 0))[0] return peak_indices def find_valleys_by_derivative(wavenumber, reflectivity, threshold=1e-3): first_derivative = np.gradient(reflectivity, wavenumber) second_derivative = np.gradient(first_derivative, wavenumber) valley_indices = np.where((np.abs(first_derivative) < threshold) & (second_derivative > 0))[0] return valley_indices # 4. 提取主频信息 def get_dominant_frequency(signal, fs=1.0): n = len(signal) fft_result = np.fft.fft(signal) freqs = np.fft.fftfreq(n, d=1/fs) magnitude = np.abs(fft_result) dominant_freq = freqs[np.argmax(magnitude[1:n//2]) + 1] return abs(dominant_freq) def plot_fft(signal, fs=1.0, title="FFT"): n = len(signal) fft_result = np.fft.fft(signal) freqs = np.fft.fftfreq(n, d=1/fs) magnitude = np.abs(fft_result) fig = plt.figure(figsize=(10, 4)) plt.plot(freqs[:n//2], magnitude[:n//2], label='幅值') plt.title(title, fontsize=14) plt.xlabel('频率(cycles per cm⁻¹)', fontsize=12) plt.ylabel('幅值', fontsize=12) plt.grid(True) plt.legend() plt.tight_layout() plt.close() return fig # 5. 数据保存函数 def save_peaks_valleys(save_path, wavenumber, reflectivity, peaks, valleys): peaks_data = pd.DataFrame({ '波数': wavenumber[peaks], '反射率': reflectivity[peaks] }) valleys_data = pd.DataFrame({ '波数': wavenumber[valleys], '反射率': reflectivity[valleys] }) peaks_data.to_csv(os.path.join(save_path, '波峰.csv'), index=False) valleys_data.to_csv(os.path.join(save_path, '波谷.csv'), index=False) def save_reflectivity_data(save_path, wavenumber, reflectivity): df = pd.DataFrame({ '波数': wavenumber, '反射率': reflectivity }) df.to_csv(os.path.join(save_path, '去噪后反射率数据.csv'), index=False) def save_dominant_frequency(save_path, freq): with open(os.path.join(save_path, '主频信息.txt'), 'w', encoding='utf-8') as f: f.write(f"信号主频(频域): {freq:.4f} cycles per cm⁻¹\n") def save_data_stats(save_path, file_name, negative_num, over_one_num): with open(os.path.join(save_path, '原始数据统计.txt'), 'w', encoding='utf-8') as f: f.write(f"文件名: {file_name}\n") f.write(f"反射率小于0的个数: {negative_num}\n") f.write(f"反射率大于100的个数: {over_one_num}\n") def save_log(save_path, log_content): with open(os.path.join(save_path, '日志.txt'), 'w', encoding='utf-8') as f: f.write(log_content) # 新增函数:波数转波长 def wavenumber_to_wavelength(wavenumber): return 1e4 / wavenumber # cm⁻¹ -> μm # 新增函数:折射率经验公式 def refractive_index(wavelength): term1 = (0.20075 * wavelength ** 2) / (wavelength ** 2 + 12.07224) term2 = (5.54861 * wavelength ** 2) / (wavelength ** 2 - 0.02641) term3 = (35.65066 * wavelength ** 2) / (wavelength ** 2 - 1268.24708) return np.sqrt(1 + term1 + term2 + term3) #复折射率 def refractive_index_complex(wavelength, reflectance): """ 计算复折射率 n = n0 + i*kappa :param wavelength: 波长 (μm) :param reflectance: 反射率 (0~100) :return: 复折射率 """ n0 = refractive_index(wavelength) # 实部 R = reflectance / 100.0 # 归一化反射率 try: kappa = np.sqrt((R * (n0 + 1) ** 2 - (n0 - 1) ** 2) / (1 - R)) except: kappa = 0.0 # 防止负数开根号 return complex(n0, kappa) # 新增函数:厚度计算 def calculate_thickness(P, wavelength, theta1_deg): theta1 = np.deg2rad(theta1_deg) sin_theta1 = np.sin(theta1) n = refractive_index(wavelength) sqrt_n2_sin2 = np.sqrt(n ** 2 - sin_theta1 ** 2) d = ((P - 0.5) * wavelength * sqrt_n2_sin2) / (2 * n ** 2) return d # 新增函数:从波峰波谷数据计算干涉周期 def compute_interference_period_weighted(peaks_wavenumber, valleys_wavenumber): """ 使用多对波峰波谷计算干涉周期 P,取加权平均 :param peaks_wavenumber: 波峰波数数组 :param valleys_wavenumber: 波谷波数数组 :return: 干涉周期 P """ if len(peaks_wavenumber) < 2 or len(valleys_wavenumber) < 2: return None peak_diffs = np.diff(peaks_wavenumber) valley_diffs = np.diff(valleys_wavenumber) all_diffs = np.concatenate([peak_diffs, valley_diffs]) # 使用中位数代替均值,提高稳健性 P = np.median(all_diffs) return P def calculate_thickness_complex(P, wavelength, theta1_deg, reflectance): """ 使用复折射率计算厚度 :param P: 干涉周期 :param wavelength: 波长 (μm) :param theta1_deg: 入射角 (度) :param reflectance: 反射率 (%) :return: 厚度 d (μm) """ theta1 = np.deg2rad(theta1_deg) sin_theta1 = np.sin(theta1) n_complex = refractive_index_complex(wavelength, reflectance) n0 = n_complex.real kappa = n_complex.imag # 使用复折射率实部计算 sqrt(n0^2 - sin^2θ1) sqrt_n2_sin2 = np.sqrt(n0 ** 2 - sin_theta1 ** 2) # 考虑吸收效应修正 absorption_correction = 1.0 + (kappa / n0) ** 2 d = ((P - 0.5) * wavelength * sqrt_n2_sin2) / (2 * n0 ** 2) * absorption_correction return d # 新增函数:计算并保存厚度 def compute_and_save_thickness(save_path, peaks_wavenumber, valleys_wavenumber, reflectivity_denoised, theta1_deg): if len(peaks_wavenumber) < 2 or len(valleys_wavenumber) < 2: print("无法计算干涉周期,波峰/波谷数量不足") return P = compute_interference_period_weighted(peaks_wavenumber, valleys_wavenumber) if P is None: print("无法提取干涉周期") return # 取波峰波谷之间的平均波长 all_wavenumbers = np.concatenate([peaks_wavenumber, valleys_wavenumber]) mean_wavenumber = np.mean(all_wavenumbers) wavelength = wavenumber_to_wavelength(mean_wavenumber) # 取平均反射率用于复折射率计算 mean_reflectance = np.mean(reflectivity_denoised) # 计算厚度(使用复折射率) thickness = calculate_thickness_complex(P, wavelength, theta1_deg, mean_reflectance) # 保存结果 result_file = os.path.join(save_path, '厚度计算结果.txt') with open(result_file, 'w', encoding='utf-8') as f: f.write(f"入射角: {theta1_deg}°\n") f.write(f"干涉周期: {P:.4f} cm⁻¹\n") f.write(f"平均波长: {wavelength:.4f} μm\n") f.write(f"平均反射率: {mean_reflectance:.2f}%\n") f.write(f"计算厚度: {thickness:.4f} μm\n") f.write(f"公式:d = (P - 0.5) * λ * sqrt(n₀² - sin²θ1) / (2n₀²) * (1 + (κ/n₀)²)\n") print(f"\n厚度计算结果已保存至:{result_file}") print(f"厚度:{thickness:.4f} μm") # 6. 主程序 if __name__ == "__main__": result_folder = create_result_folder("results") file_paths = [ os.path.join(path, "附件1.xlsx"), os.path.join(path, "附件2.xlsx") ] error_report = {} for file in file_paths: log_output = io.StringIO() console = sys.stdout sys.stdout = log_output print(f"\n对于文件:{os.path.basename(file)}") wavenumber, reflectivity = get_data(file) if wavenumber is not None and reflectivity is not None: file_name = os.path.splitext(os.path.basename(file))[0] save_path = os.path.join(result_folder, file_name) create_result_folder(save_path) # 保存原始数据统计 file_name_full = os.path.basename(file) negative_num = np.sum(reflectivity < 0) over_one_num = np.sum(reflectivity > 100) save_data_stats(save_path, file_name_full, negative_num, over_one_num) # 数据预处理 wavenumber_clean, reflectivity_clean = outlier_removal(wavenumber, reflectivity) reflectivity_smooth = apply_savgol_filter(reflectivity_clean) reflectivity_denoised = wavelet_denoise(reflectivity_smooth, mode_ext='constant') reflectivity_denoised = equal_length(reflectivity_denoised, len(wavenumber_clean)) reflectivity_normalized = normalize(reflectivity_denoised) # 去噪前后对比图 plt.figure(figsize=(12, 6)) plt.plot(wavenumber_clean, reflectivity_clean, label='原始数据', alpha=0.6, color='gray') plt.plot(wavenumber_clean, reflectivity_denoised, label='去噪后数据', color='blue') plt.title(f'{file_name} 去噪前后对比', fontsize=14) plt.xlabel('波数 (cm⁻¹)') plt.ylabel('反射率 (%)') plt.legend() plt.grid(True) plt.gca().invert_xaxis() plt.tight_layout() save_plot(save_path, '去噪前后对比图.png') # 波峰波谷检测(导数法) peak_indices = find_peaks_by_derivative(wavenumber_clean, reflectivity_denoised, threshold=1e-3) valley_indices = find_valleys_by_derivative(wavenumber_clean, reflectivity_denoised, threshold=1e-3) plt.figure(figsize=(12, 6)) plt.plot(wavenumber_clean, reflectivity_denoised, label='反射率', color='blue') plt.scatter(wavenumber_clean[peak_indices], reflectivity_denoised[peak_indices], marker='o', s=30, color='red', label='波峰(极大值)') plt.scatter(wavenumber_clean[valley_indices], reflectivity_denoised[valley_indices], marker='x', s=30, color='green', label='波谷(极小值)') plt.title('导数法检测干涉极大值点与极小值点', fontsize=14) plt.xlabel('波数 (cm⁻¹)', fontsize=12) plt.ylabel('反射率 (%)', fontsize=12) plt.legend() plt.grid(True) plt.gca().invert_xaxis() plt.tight_layout() save_plot(save_path, '导数法波峰波谷检测图.png') save_peaks_valleys(save_path, wavenumber_clean, reflectivity_denoised, peak_indices, valley_indices) save_reflectivity_data(save_path, wavenumber_clean, reflectivity_denoised) # 在主程序中,找到波峰波谷后插入以下代码: peak_wavenumbers = wavenumber_clean[peak_indices] valley_wavenumbers = wavenumber_clean[valley_indices] # 获取入射角(根据文件名判断) if "附件1" in file_name: theta1_deg = 10 elif "附件2" in file_name: theta1_deg = 15 else: theta1_deg = 0 # 调用厚度计算函数 compute_and_save_thickness(save_path, peak_wavenumbers, valley_wavenumbers, reflectivity_denoised, theta1_deg) # 主频信息 delta_wavenumber = np.mean(np.diff(wavenumber_clean)) fs = 1.0 / delta_wavenumber dominant_freq = get_dominant_frequency(reflectivity_denoised, fs=fs) save_dominant_frequency(save_path, dominant_freq) # 频域图 fig = plot_fft(reflectivity_clean, fs=fs, title='原始信号频域分析') save_plot(save_path, '原始信号频域图.png', fig) fig = plot_fft(reflectivity_denoised, fs=fs, title='去噪后信号频域分析') save_plot(save_path, '去噪后信号频域图.png', fig) # 误差分析 plt.figure(figsize=(12, 6)) plt.plot(wavenumber_clean, reflectivity_clean, label='原始信号', alpha=0.6, color='gray') plt.plot(wavenumber_clean, reflectivity_denoised, label='去噪后信号', color='blue') plt.plot(wavenumber_clean, reflectivity_clean - reflectivity_denoised, label='残差', linestyle='--', color='red') plt.title('去噪前后误差分析') plt.xlabel('波数 (cm⁻¹)') plt.ylabel('反射率 (%)') plt.legend() plt.grid(True) plt.gca().invert_xaxis() plt.tight_layout() save_plot(save_path, '去噪前后误差分析图.png') print(f"\n【主频稳定性分析 - {file_name}】") freqs = [] for _ in range(10): noisy_signal = reflectivity_denoised + np.random.normal(0, 0.1, len(reflectivity_denoised)) freq = get_dominant_frequency(noisy_signal, fs=fs) freqs.append(freq) mean_freq = np.mean(freqs) std_freq = np.std(freqs) print(f"主频估计均值:{mean_freq:.6f}±{std_freq:.6f}") plt.hist(freqs, bins=10, edgecolor='black') plt.title('主频估计稳定性分析') plt.xlabel('主频 (cycles per cm⁻¹)') plt.ylabel('频次') plt.grid(True) save_plot(save_path, '主频估计稳定性分析图.png') # 收集误差数据 error_report[file_name] = { "去噪前后信号MSE": np.mean((reflectivity_clean - reflectivity_denoised) ** 2), "主频估计均值": float(mean_freq), "主频估计标准差": float(std_freq), "波峰检测数量": int(len(peak_indices)), "波谷检测数量": int(len(valley_indices)) } # 恢复标准输出并保存日志 sys.stdout = console save_log(save_path, log_output.getvalue()) # 保存误差分析报告 with open(os.path.join(result_folder, "误差分析报告.json"), "w", encoding="utf-8") as f: json.dump(error_report, f, ensure_ascii=False, indent=4) print("误差分析报告已保存至:", os.path.join(result_folder, "误差分析报告.json"))
最新发布
09-06
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值