一、插值

         为了降低实时处理时间，采用粗暴地边读边插值的方法，如果当前数据包间隔上一个数据包超过100，就插值。

for i = 1:row
    %读取第i条数据
    csi_entry_N = csi_trace{i};%结构体形式
    if isstruct(csi_entry_N)  == 0 %判断是否为空，0为空
       continue;
    end

    %读取并存储第i条数据的时间戳
    current_timestamp = csi_trace{i}.timestamp_low;
    result_matrix(i,241) = csi_trace{i}.timestamp_low;
    
    % 检测时间戳间隔
    if i > 1 % 添加这个条件来避免第一次循环时的错误
        timestamp_interval = current_timestamp - prev_timestamp;
 
        % 计算需要补足的数据包数量
        num_missing_packets = round(timestamp_interval / 100) - 1;
 
        % 使用矩阵运算进行线性插值
        if num_missing_packets > 0
            prev_amplitude = result_matrix(i-1, 1:120);
            prev_phase = result_matrix(i-1, 121:240);
            current_amplitude = result_matrix(i, 1:120);
            current_phase = result_matrix(i, 121:240);
 
            % 计算振幅和相位的差值
            amplitude_diff = current_amplitude - prev_amplitude;
            phase_diff = current_phase - prev_phase;
 
            % 计算插值的倍数
            interpolation_factors = (1:num_missing_packets) / (num_missing_packets + 1);
 
            % 广播插值的倍数，并计算插值后的振幅和相位
            interpolated_amplitude = prev_amplitude + amplitude_diff .* interpolation_factors';
            interpolated_phase = prev_phase + phase_diff .* interpolation_factors';
 
            % 插入补足的数据包
            new_row_index_range = new_row_index + (1:num_missing_packets);
            result_matrix(new_row_index_range, 1:120) = interpolated_amplitude;
            result_matrix(new_row_index_range, 121:240) = interpolated_phase;
            result_matrix(new_row_index_range, 241) = prev_timestamp + (1:num_missing_packets) * 100; % 更新时间戳
 
            % 更新新数据包的索引
            new_row_index = new_row_index + num_missing_packets;
        end
    end
     % 更新前一个时间戳
    prev_timestamp = current_timestamp;

    %读取第i条数据的CSI值
    csi = get_scaled_csi(csi_entry_N);%复数矩阵,2*2*30
    csi1 = squeeze(csi);% 去掉全为1的数组，复数矩阵
    [h,l,y]=size(csi);H=[H h];L =[L l];
    if (h==2 && l==2)
        % 提取当前时刻的信号数据包
        current_packet = csi1;
 
        % 计算振幅信息
        amplitude = abs(current_packet);
        %disp(amplitude);
        % 计算相位信息
        phase = angle(current_packet);
        
        %disp("------------------------")
        % 将振幅和相位信息存储到result_matrix中
        result_matrix(i, 1:120) = reshape(amplitude, 1, []); %把amplitude这个数组改编为1*n的形式
        %disp(result_matrix);
        result_matrix(i, 121:240) = reshape(phase, 1, []);

        % 更新新数据包的索引
        new_row_index = new_row_index + 1;
    else
        continue;
    end
    

end

二、hanple极值点去除

% 这一节得到是降噪后的result_matrix
% 需要设置的参数：window_size、num_dev
% 获取振幅数据
amplitude_data = result_matrix(:, 1:120);
 
% 可视化
figure;
% 画出原始数据
subplot(2, 1, 1);  % 2x1的子图布局，选中第1个子图
plot(amplitude_data);
xlabel('时间');
ylabel('振幅');
title('振幅随时间的变化（原始）');
hold on
 
% 设定hampel滤波器参数
window_size = 5; % 窗口大小，可根据需要调整
num_dev = 1; % 定义异常值的标准偏差范围
 
% 对每一个子载波进行处理
for k = 1:90
    % 使用hampel函数进行滤波
    filtered_data = hampel(amplitude_data(:, k), window_size, num_dev);
    
    % 存回到result_matrix
    result_matrix(:, k) = filtered_data;
end
 
% 画出经Hampel处理后的数据
subplot(2, 1, 2);  % 2x1的子图布局，选中第2个子图
plot(result_matrix(:, 1:120));
xlabel('时间');
ylabel('振幅');
title('振幅随时间的变化（Hampel）');

三、低通滤波器处理

%% 降噪二：低通滤波器处理
% 这一节得到是降噪后的result_matrix
% 需要设置的参数：Fc、order、Fs
 
% 提取振幅和相位数据
amplitude_data = result_matrix(:, 1:120);
 
% 可视化
figure;
% 画出原始数据
subplot(2, 1, 1);  % 2x1的子图布局，选中第1个子图
plot(amplitude_data);
xlabel('时间');
ylabel('振幅');
title('振幅随时间的变化（原始）');
hold on
 
% 设计低通滤波器
Fc = 40; % 截止频率，根据实际情况进行调整，但确保小于15Hz
order = 4; % 滤波器阶数，可以根据需要进行调整
Fs = 100; % 采样频率，已经给出是30Hz
lowpass_filter = designfilt('lowpassiir', 'FilterOrder', order, 'HalfPowerFrequency', Fc, 'SampleRate', Fs, 'DesignMethod', 'butter');
 
% 对每个子载波上的振幅和相位数据应用低通滤波器
filtered_data = zeros(size(amplitude_data));
for i = 1:120
    filtered_data(:, i) = filtfilt(lowpass_filter, amplitude_data(:, i));
end
 
% 将滤波后的振幅和相位数据保存回result_matrix
result_matrix(:, 1:120) = filtered_data;
 
% 画出经低通滤波器处理后的数据
subplot(2, 1, 2);  % 2x1的子图布局，选中第2个子图
plot(result_matrix(:, 1:120));
xlabel('时间');
ylabel('振幅');
title('振幅随时间的变化（lowpass）');
hold off

四、离散小波变换

%% 降噪三：离散小波变换
% 这一节得到是降噪后的result_matrix
% 提取振幅和相位数据
amplitude_data = result_matrix(:, 1:120);
% 可视化
figure;
% 画出原始数据
subplot(2, 1, 1);  % 2x1的子图布局，选中第1个子图
plot(amplitude_data);
xlabel('时间');
ylabel('振幅');
title('振幅随时间的变化（原始）');
hold on

% 初始化去噪后的数据矩阵
filtered_data = zeros(size(amplitude_data));

% 对每个子载波上的振幅和相位数据应用小波去噪
for i = 1:120
    maxLevel = wmaxlev(size(amplitude_data, 1), 'sym4'); % 获取最大可能的分解级别
    level = min(10, maxLevel); % 保证选择的分解级别不超过最大可能级别
    filtered_data(:, i) = wden(amplitude_data(:, i), 'heursure', 's', 'one', level, 'sym4');
end

% 将去噪后的振幅和相位数据保存回result_matrix
result_matrix(:, 1:120) = filtered_data;

% 画出经小波变换处理后的数据
subplot(2, 1, 2);  % 2x1的子图布局，选中第2个子图
plot(result_matrix(:, 1:120));
xlabel('时间');
ylabel('振幅');
title('振幅随时间的变化（DWT）');
hold off