%% 1. 数据读取与预处理
full_data = readmatrix('dingji.csv');
time1 = full_data(:, 1);
speed1 = -full_data(:, 2) / 180 * pi;
% 可视化
figure;
plot(time1, speed1, 'r-', 'LineWidth', 1.2);
xlabel('时间 (秒)');
ylabel('速度');
title('原始速度数据');
grid on;
%% 1. 数据生成替代 CSV 读取
% 参数设置
lb = [0.680, 1.784, -pi, 0.945, -inf]; % 下界
ub = [1.145, 2.15, pi, 2.09 - 0.68, inf]; % 上界
% 固定参数或随机选择：
A = 0.82; % 振幅，在 lb(1) 和 ub(1) 之间
omega = 1.92; % 频率，在 lb(2) 和 ub(2) 之间
phi = 0.1; % 相位，在 lb(3) 和 ub(3) 之间
v0 = 1.27; % 初始速度，在 lb(4) 和 ub(4) 之间
s0 = 0.2; % 初始位移
% 时间设置
num_points = 1001;
dt = 0.01; % 采样频率为 0.01 秒
time = (0:num_points-1)' * dt;
% 生成速度数据
speed = -A ./ omega .* cos(omega .* time + phi) + v0 .* time + s0;
% 可选：添加一点噪声以更贴近真实测量数据
noise_level = 0.03; % 噪声幅度
speed = speed + noise_level * randn(size(speed));
% 显示前几个数据点验证
disp('前10个时间点:');
disp(time(1:10));
disp('前10个速度点:');
disp(speed(1:10));
% 可视化
figure;
plot(time, speed, 'b-', 'LineWidth', 1.2);
xlabel('时间 (秒)');
ylabel('速度');
title('合成生成的速度数据');
grid on;
%% 2. 计算加速度（使用中心差分）
% 原始数据 speed1 的加速度
dt1 = diff(time1); % 各个时间间隔
valid_idx1 = (dt1(1:end-1) + dt1(2:end)) ~= 0; % 避免除以0
acc1 = zeros(length(speed1) - 2, 1);
time_acc1 = zeros(length(speed1) - 2, 1);
for i = 2:length(speed1)-1
delta_t = time1(i+1) - time1(i-1);
if delta_t == 0
acc1(i-1) = NaN;
else
acc1(i-1) = (speed1(i+1) - speed1(i-1)) / delta_t;
end
time_acc1(i-1) = (time1(i+1) + time1(i-1)) / 2;
end
% 合成数据 speed 的加速度
dt2 = diff(time); % 各个时间间隔
valid_idx2 = (dt2(1:end-1) + dt2(2:end)) ~= 0;
acc2 = zeros(length(speed) - 2, 1);
time_acc2 = zeros(length(speed) - 2, 1);
for i = 2:length(speed)-1
delta_t = time(i+1) - time(i-1);
if delta_t == 0
acc2(i-1) = NaN;
else
acc2(i-1) = (speed(i+1) - speed(i-1)) / delta_t;
end
time_acc2(i-1) = (time(i+1) + time(i-1)) / 2;
end
%% 3. 可视化加速度对比图
figure;
subplot(2, 1, 1);
plot(time_acc1, acc1, 'r-', 'LineWidth', 1.2);
xlabel('时间 (秒)');
ylabel('加速度');
title('原始数据计算出的加速度');
grid on;
subplot(2, 1, 2);
plot(time_acc2, acc2, 'b-', 'LineWidth', 1.2);
xlabel('时间 (秒)');
ylabel('加速度');
title('合成数据计算出的加速度');
grid on;
%% 4. 对加速度进行高斯滤波
% 设置高斯滤波参数
sigma =20; % 高斯核标准差，越大越平滑
window_size = ceil(6*sigma); % 初始窗口大小
if mod(window_size, 2) == 0 % 确保是奇数
window_size = window_size + 1;
end
window_size = max(3, window_size); % 确保至少有3个点
% 对原始数据的加速度进行高斯滤波
acc1_smooth = imgaussfilt(acc1, sigma, 'FilterSize', window_size);
% 对合成数据的加速度进行高斯滤波
acc2_smooth = imgaussfilt(acc2, sigma, 'FilterSize', window_size);
%% 5. 可视化对比图
figure;
% 原始数据加速度对比
subplot(2, 1, 1);
plot(time_acc1, acc1, 'r-', 'LineWidth', 1.0);
hold on;
plot(time_acc1, acc1_smooth, 'k-', 'LineWidth', 1.5);
legend('原始加速度', '高斯滤波结果');
title('原始数据加速度与滤波对比');
xlabel('时间 (秒)');
ylabel('加速度');
grid on;
hold off;
% 合成数据加速度对比
subplot(2, 1, 2);
plot(time_acc2, acc2, 'b-', 'LineWidth', 1.0);
hold on;
plot(time_acc2, acc2_smooth, 'k-', 'LineWidth', 1.5);
legend('原始加速度', '高斯滤波结果');
title('合成数据加速度与滤波对比');
xlabel('时间 (秒)');
ylabel('加速度');
grid on;
hold off;
请对滤波后的加速度曲线进行拟合， % 初始参数估计（粗略估计频率和振幅）
initial_a = (0.780 + 1.045) / 2;
initial_w = (1.884 + 2) / 2; % 基于采样时间范围估计频率
initial_w0 = 0;
initial_b = 2.09 - initial_a;
initial_params = [initial_a, initial_w, initial_w0, initial_b];
% 设置参数下界：a必须为正，其他参数无限制
lb = [0.780, 1.045, -pi]; % a的下界设为1e-6避免为0
ub = [1.045, 2, pi]; % 所有参数无上界
% 使用最小二乘法拟合正弦函数
try
% 定义拟合函数
sine_model = @(params, t) params(1) * sin(params(2) * t + params(3)) + 2.09 - params(1);
% 非线性最小二乘拟合
fitted_params = lsqcurvefit(sine_model, initial_params, t_sample, y_sample, lb, ub);
catch
continue; % 拟合失败则跳过当前迭代
end
最新发布
扫一扫
1100
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
