用于分析介电谱数据的MATLAB项目

文件结构图

plaintext

├── main.m                     % 主脚本：执行全流程分析
├── functions/
│   ├── generate_simulated_data.m % 生成模拟介电谱数据
│   ├── fit_spectrum.m          % 多项式拟合函数
│   ├── calculate_slope_function.m % 计算斜率函数
│   ├── find_characteristic_frequency.m % 寻找特征频率点
│   ├── calculate_activation_energy.m % 计算活化能
│   ├── normalize_spectrum.m      % 频谱归算函数
├── plots/                      % 生成的图片存储目录
│   ├── dielectric_loss_curves.png
│   ├──拟合曲线.png
│   ├──特征频率点.png
│   ├──活化能与含水率关系.png
│   ├──归算曲线对比.png
│   └──复介电常数归算.png
└── parameters.m                % 模型参数定义

1. 参数定义文件（parameters.m）

% 定义实验参数
params.temperatures = [-40, -20, 10]; % 修改后的温度范围（℃）
params.moistures = [0.2, 1.2, 2.7];   % 含水率（%）
params.reference_temp = 5;            % 参考温度（℃，用于验证）
params.frequency_range = logspace(-3, 4, 100); % 频率范围（Hz）
params.poly_order = 4;                % 多项式拟合阶数
params.boltzmann = 1.38e-23;          % 玻尔兹曼常数（J/K）

2. 模拟数据生成函数（functions/generate_simulated_data.m）

function [tan_delta, epsilon_real, epsilon_imag] = generate_simulated_data(freq, temp, moisture)
% 基于论文模型生成模拟介电损耗和复介电常数
% 注：实际应用需替换为真实实验数据

% 弛豫时间活化能与含水率的指数关系（论文式11）
E_tau = 0.00841 * exp(moisture / 0.00964) + 0.70956; % eV，转换为J需乘以1.6e-19
E_tau_J = E_tau * 1.6e-19;

% 弛豫时间与温度的关系（论文式7）
tau = @(T) exp(E_tau_J ./ (params.boltzmann * (T + 273.15))); % 假设tau0=1s

% 模拟介电损耗因数（简化Debye模型）
tan_delta = zeros(size(freq));
epsilon_real = zeros(size(freq));
epsilon_imag = zeros(size(freq));

for i = 1:length(freq)
    tau_val = tau(temp);
    omega = 2*pi*freq(i);
    epsilon_s = 3 + 0.5*moisture; % 静态介电常数（模拟值）
    epsilon_inf = 2;               % 光频介电常数（模拟值）
    epsilon_real(i) = epsilon_inf + (epsilon_s - epsilon_inf) ./ (1 + (omega*tau_val).^2);
    epsilon_imag(i) = (epsilon_s - epsilon_inf)*omega*tau_val ./ (1 + (omega*tau_val).^2);
    tan_delta(i) = epsilon_imag(i) ./ epsilon_real(i);
end

3. 主脚本（main.m）

%% 清空环境并加载参数
clear; clc;
run('parameters.m');
mkdir('plots'); % 创建图片存储目录

%% 生成模拟数据并绘图
figure('Position', [100 100 800 600]);
for m = 1:length(params.moistures)
    moisture = params.moistures(m);
    for t = 1:length(params.temperatures)
        temp = params.temperatures(t);
        [tan_delta, epsilon_real, epsilon_imag] = functions.generate_simulated_data(params.frequency_range, temp, moisture);
        
        % 绘制介电损耗曲线
        subplot(2,2,m);
        semilogx(params.frequency_range, tan_delta, 'o-', 'DisplayName', sprintf('T=%d℃', temp));
        title(sprintf('含水率%.1f%%', moisture));
        xlabel('频率 (Hz)'); ylabel('tanδ');
        legend('Location', 'best');
    end
end
saveas(gcf, 'plots/dielectric_loss_curves.png'); % 保存图1：不同含水率变温介电损耗曲线


%% 多项式拟合与特征频率提取
for m = 1:length(params.moistures)
    moisture = params.moistures(m);
    for t = 1:length(params.temperatures)
        temp = params.temperatures(t);
        [tan_delta, ~, ~] = functions.generate_simulated_data(params.frequency_range, temp, moisture);
        
        % 对数化处理（论文式6前处理）
        log_freq = log10(params.frequency_range);
        log_tan_delta = log10(tan_delta);
        
        % 多项式拟合
        coeff = polyfit(log_freq, log_tan_delta, params.poly_order);
        fit_tan_delta = 10.^polyval(coeff, log_freq);
        
        % 绘制拟合曲线
        figure;
        semilogx(params.frequency_range, tan_delta, 'o', params.frequency_range, fit_tan_delta, '-');
        title(sprintf('含水率%.1f%%, T=%d℃ 拟合曲线', moisture, temp));
        xlabel('频率 (Hz)'); ylabel('tanδ');
        saveas(gcf, sprintf('plots/拟合曲线_%.1f%%_T%d℃.png', moisture, temp)); % 保存拟合曲线（3-4个图）
    end
end


%% 计算特征频率点与活化能
activation_energy = zeros(length(params.moistures), 1);
for m = 1:length(params.moistures)
    moisture = params.moistures(m);
    freqs = zeros(length(params.temperatures), 1);
    
    for t = 1:length(params.temperatures)
        temp = params.temperatures(t);
        [tan_delta, ~, ~] = functions.generate_simulated_data(params.frequency_range, temp, moisture);
        log_freq = log10(params.frequency_range);
        log_tan_delta = log10(tan_delta);
        
        % 计算斜率函数（论文式6）
        slope = functions.calculate_slope_function(log_freq, log_tan_delta);
        
        % 寻找极小值点（特征频率）
        [~, min_idx] = min(slope);
        freqs(t) = params.frequency_range(min_idx);
        
        % 绘制斜率函数与特征点
        figure;
        plot(log_freq, slope);
        hold on;
        plot(log_freq(min_idx), slope(min_idx), 'ro', 'MarkerSize', 10);
        title(sprintf('含水率%.1f%%, T=%d℃ 特征频率点', moisture, temp));
        xlabel('log10(频率)'); ylabel('d(log10(tanδ))/d(log10(f))');
        saveas(gcf, sprintf('plots/特征频率点_%.1f%%_T%d℃.png', moisture, temp)); % 保存特征频率图（3个图）
    end
    
    % 计算活化能（论文式9）
    T1 = params.temperatures(1) + 273.15; % K
    T2 = params.temperatures(2) + 273.15;
    f1 = log(freqs(1));
    f2 = log(freqs(2));
    activation_energy(m) = params.boltzmann * (f2 - f1) / (1/T1 - 1/T2) / 1.6e-19; % 转换为eV
end


%% 绘制活化能与含水率关系
figure;
moistures = params.moistures;
plot(moistures, activation_energy, 'o-', 'LineWidth', 2);
xlabel('含水率 (%)'); ylabel('弛豫时间活化能 (eV)');
title('活化能与含水率关系');
% 拟合指数曲线（论文式11）
fit = fit(moistures', activation_energy', 'exp1');
hold on;
fplot(@(x) fit.a*exp(x/fit.b)+fit.c, [0, 3.5], 'r--');
legend('实测数据', '拟合曲线');
saveas(gcf, 'plots/活化能与含水率关系.png'); % 保存图5


%% 频谱归算至参考温度（5℃）
reference_T = params.reference_temp + 273.15; % K
for m = 1:length(params.moistures)
    moisture = params.moistures(m);
    E_tau = 0.00841 * exp(moisture / 0.00964) + 0.70956; % 论文式11
    E_tau_J = E_tau * 1.6e-19;
    
    figure;
    for t = 1:length(params.temperatures)
        temp = params.temperatures(t);
        T_kelvin = temp + 273.15;
        [tan_delta, epsilon_real, epsilon_imag] = functions.generate_simulated_data(params.frequency_range, temp, moisture);
        
        % 归算频率（论文式12）
        normalized_freq = params.frequency_range .* exp( E_tau_J/params.boltzmann * (1/reference_T - 1/T_kelvin) );
        
        % 绘制归算后曲线
        semilogx(normalized_freq, tan_delta, '.-', 'DisplayName', sprintf('归算自%d℃', temp));
    end
    
    % 绘制参考温度下的"实测数据"（模拟5℃数据）
    [tan_delta_ref, ~, ~] = functions.generate_simulated_data(params.frequency_range, params.reference_temp, moisture);
    semilogx(params.frequency_range, tan_delta_ref, 'k--', 'LineWidth', 2, 'DisplayName', '5℃实测数据');
    
    title(sprintf('含水率%.1f%% 归算至5℃', moisture));
    xlabel('频率 (Hz)'); ylabel('tanδ');
    legend('Location', 'best');
    saveas(gcf, sprintf('plots/归算曲线对比_%.1f%%.png', moisture)); % 保存归算验证图（3个图）
end

4. 辅助函数（functions / 目录下）

4.1 计算斜率函数（calculate_slope_function.m）

function slope = calculate_slope_function(log_freq, log_tan_delta)
% 计算对数坐标系下的斜率（论文式6）
slope = diff(log_tan_delta) ./ diff(log_freq);
% 插值到原始频率点
slope = [slope(1); (slope(1:end-1)+slope(2:end))/2; slope(end)];
end

4.2 寻找特征频率点（find_characteristic_frequency.m）

function [f_char] = find_characteristic_frequency(log_freq, slope)
% 寻找斜率函数的极小值点作为特征频率
[~, min_idx] = min(slope);
f_char = 10^log_freq(min_idx);
end

4.3 归算函数（normalize_spectrum.m）

function [normalized_freq, normalized_data] = normalize_spectrum(freq, data, temp, moisture, reference_temp)
% 频谱归算核心函数（论文式12）
E_tau = 0.00841 * exp(moisture / 0.00964) + 0.70956; % eV
E_tau_J = E_tau * 1.6e-19;
T = temp + 273.15;
T_ref = reference_temp + 273.15;

normalized_freq = freq .* exp( E_tau_J/params.boltzmann * (1/T_ref - 1/T) );
normalized_data = interp1(freq, data, normalized_freq, 'spline');
end

生成的图片示例（6-7 个图）

1. 图 1：不同含水率变温介电损耗曲线
   （3 个子图，每个含水率对应 - 40℃、-20℃、10℃曲线）

2. 拟合曲线
   （每个含水率 + 温度组合 1 个图，共 3×3=9 个，可按需选择 3 个展示）

3. 特征频率点
   （每个含水率 + 温度组合 1 个图，共 3×3=9 个，可按需选择 3 个展示）

4. 图 5：活化能与含水率关系
   （含指数拟合曲线）

5. 归算曲线对比
   （每个含水率 1 个图，展示归算至 5℃后与模拟 5℃数据的重合度，共 3 个图）

注意事项

1. 实际数据替换：请将generate_simulated_data函数替换为真实实验数据读取代码。
2. 低温验证：论文中未涉及低温场景，需注意低温下油纸绝缘可能出现的物理变化（如油凝固），建议先通过实验验证模型适用性。
3. 参数调整：多项式拟合阶数、频率范围等可根据实际数据调整。

 

运行主脚本main.m即可生成所有图形和数据，结果存储在plots/目录下。