6-27 gauge数据处理

clear
clc
close all

%% Case parameters
T = 1.4; % [s]
A = 0.145; % [m]
d = 1.25; % m

g = 9.81; % m/s^2

% Sample frequencies - Radar
SR = 10;

% Sample frequency - Gauge
SR_gauge = 128;

% Types of roughness
roughness_type = {'Fan'}; 

% Number of roughness types
nType = length(roughness_type);

%% Adding paths to any relevant folders
% Generic MATLAB functions
addpath('/Users/liuxinyu/Library/CloudStorage/OneDrive-ImperialCollegeLondon/thesis/scripts/WaveGenTutorial/Analysis/2.Functions');
% Spiking toolbox
addpath('/Users/liuxinyu/Library/CloudStorage/OneDrive-ImperialCollegeLondon/thesis/scripts/Analysis files/despiking_toolbox');
% Non-linear Regular Wave functions
addpath('/Users/liuxinyu/Library/CloudStorage/OneDrive-ImperialCollegeLondon/thesis/scripts/Analysis files/Non-linear Regular');


%% READ GAUGE DATA
folders_gauges = {'/Users/liuxinyu/Library/CloudStorage/OneDrive-ImperialCollegeLondon/thesis/scripts/data_senthuren/Data/Middle Gauge/RegularWaves/Output/mat',
    '/Users/liuxinyu/Library/CloudStorage/OneDrive-ImperialCollegeLondon/thesis/scripts/data_senthuren/Data/MIROS/Roughness test/Fan/Gauge/RegularWaves/Output/mat'}

Gauges_Regular = {};
for i = 1:length(folders_gauges)
    mat_files = dir(fullfile(folders_gauges{i}, sprintf('*T%.1f_A%.3f*.mat', T, A)));

    for j = 1:length(mat_files)
        mat_file_name = mat_files(j).name;
        full_file_path = fullfile(folders_gauges{i}, mat_file_name);

        loaded_data = load(full_file_path);

        Gauges_Regular{end + 1, 1} = loaded_data;

        eta = loaded_data.eta;
        rows = (1:size(eta, 1))';

        % Interpolate NaNs per column using 'pchip'
        for col = 1:size(eta, 2)
            y = eta(:, col);
            valid_idx = ~isnan(y);
            missing_idx = isnan(y);

            % Only interpolate if there are at least 2 valid points
            if sum(valid_idx) >= 2 && any(missing_idx)
                y_interp = interp1(rows(valid_idx), y(valid_idx), rows(missing_idx), 'pchip');
                y(missing_idx) = y_interp;
                eta(:, col) = y;  % Update column
            end
        end

        Gauges_Regular{i}.eta = eta;
        gauge_eta{i} = eta;
    end
end

gauge_t = (0 : size(gauge_eta{1},1) - 1) / SR_gauge

%% 低通滤波，消除高频噪声（尤其是>20Hz的噪声）
for iCase_g = 1:length(gauge_eta)
    eta = gauge_eta{iCase_g};
    rows = (1:size(eta, 1))';

    % NaN插值
    for col = 1:size(eta, 2)
        y = eta(:, col);
        valid_idx = ~isnan(y);
        missing_idx = isnan(y);

        if sum(valid_idx) >= 2 && any(missing_idx)
            y_interp = interp1(rows(valid_idx), y(valid_idx), rows(missing_idx), 'pchip');
            y(missing_idx) = y_interp;
            eta(:, col) = y;
        end
    end

    % 更新插值后的数据
    gauge_eta{iCase_g} = eta;

    % 滤波
    f_gauge_osc = 20; % Hz
    [b,a] = butter(10, f_gauge_osc / (SR_gauge/2));

    eta_gauges_filt{iCase_g} = zeros(size(eta));
    for iGauge = 1:size(eta, 2)
        eta_gauges_filt{iCase_g}(:, iGauge) = filtfilt(b, a, eta(:, iGauge));
    end
end

%% 抗混叠滤波（Anti-Aliasing Filter）并下采样
% 抗混叠滤波设置：
% 截止频率 = 0.9 × Nyquist频率 = 0.9 × 5 Hz = 4.5 Hz

% 设置参数
SR_gauge = 128; % Gauge采样率
SR_radar = 10;  % Radar采样率

% 抗混叠滤波器设置
f_cutoff = 0.9 * (SR_radar/2); % 截止频率 = 4.5Hz
[b_aa, a_aa] = butter(10, f_cutoff / (SR_gauge/2));

% 初始化下采样数据
eta_gauges_ds = cell(1, length(gauge_eta)); 

for iCase_g = 1:length(gauge_eta)
    eta = gauge_eta{iCase_g};

    % 构造Gauge原始时间轴
    t_gauge = (0:size(eta,1)-1) / SR_gauge;
    % 构造目标时间轴（下采样后，Radar采样率）
    t_radar = (0:1/SR_radar:t_gauge(end))';

    % 抗混叠滤波
    eta_aa = zeros(size(eta));
    for iGauge = 1:size(eta,2)
        eta_aa(:, iGauge) = filtfilt(b_aa, a_aa, eta(:, iGauge));
    end

    % 下采样：使用pchip插值到t_radar
    eta_ds = zeros(length(t_radar), size(eta,2));
    for iGauge = 1:size(eta, 2)
        eta_ds(:, iGauge) = interp1(t_gauge, eta_aa(:,iGauge), t_radar, 'pchip');
    end

    % 保存
    eta_gauges_ds{iCase_g} = eta_ds;
end

%% 去均值

for iCase_g = 1:length(eta_gauges_ds)
    eta_ds = eta_gauges_ds{iCase_g};

    for iGauge = 1:size(eta_ds, 2)
        eta_gauges_ds{iCase_g}(:,iGauge) = eta_ds(:,iGauge) - mean(eta_ds(:,iGauge));
    end
end



%% analysis

%% Plot η-t for Gauge 18 in the first dataset

% 选择Gauge编号
gauge_index = 18; 

% 获取时间轴（与Radar一致）
t_gauge_ds = (0:size(eta_gauges_ds{1},1)-1) / SR; 

% 绘图
figure;
set(gcf, 'Color', 'w');
plot(t_gauge_ds, eta_gauges_ds{1}(:,gauge_index), 'b-', 'LineWidth', 1.8);

xlabel('Time [s]');
ylabel('\eta(t) [m]');
title('Wave Elevation vs Time - Gauge 18 (Set 1)');
grid on;
grid minor;
xlim([t_gauge_ds(1), t_gauge_ds(end)]);
legend('Gauge 18 - Set 1');

%%
%% Compare Gauge 18 in Set 1 and Set 2

figure;
set(gcf, 'Color', 'w');
hold on;
grid on;

plot(t_gauge_ds, eta_gauges_ds{1}(:,18), 'b-', 'LineWidth', 1.5, 'DisplayName', 'Gauge 18 - Set 1');
plot(t_gauge_ds, eta_gauges_ds{2}(:,18), 'r--', 'LineWidth', 1.5, 'DisplayName', 'Gauge 18 - Set 2');

xlabel('Time [s]');
ylabel('\eta(t) [m]');
title('Wave Elevation vs Time - Gauge 18 (Set 1 vs Set 2)');
legend('show');
xlim([t_gauge_ds(1), t_gauge_ds(end)]);
%%
%% PSD Plot for Gauge 18 (Set 1)
% linearlinear coordinates
% 数据
eta_signal = eta_gauges_ds{1}(:,18); % Set 1 中 Gauge 18

% 参数
window = hamming(256);
nfft = 1024;
noverlap = 128;
fs = SR; % 采样频率，与Radar一致（10 Hz）

% 计算PSD
[psd_gauge, f_gauge] = pwelch(eta_signal, window, noverlap, nfft, fs);

% 限制频率到0-4.5 Hz
freq_limit = 4.5;
valid_idx = f_gauge <= freq_limit;

% 绘制PSD
figure;
set(gcf, 'Color', 'w');
plot(f_gauge(valid_idx), psd_gauge(valid_idx), 'b-', 'LineWidth', 1.8);

xlabel('Frequency [Hz]');
ylabel('PSD [m^2/Hz]');
title('Power Spectral Density - Gauge 18 (Set 1)');
grid on;
grid minor;
xlim([0 freq_limit]);
ylim([0 max(psd_gauge(valid_idx))*1.1]);
legend('Gauge 18');

%%
%% PSD Plot (Log scale) for Gauge 18 (Set 1)

% 数据
eta_signal = eta_gauges_ds{1}(:,18); % Set 1 中 Gauge 18

% 参数
window = hamming(256);
nfft = 1024;
noverlap = 128;
fs = SR; % 采样频率，与Radar一致（10 Hz）

% 计算PSD
[psd_gauge, f_gauge] = pwelch(eta_signal, window, noverlap, nfft, fs);

% 限制频率到0-4.5 Hz
freq_limit = 4.5;
valid_idx = f_gauge <= freq_limit;

% 绘制PSD（对数y轴）
figure;
set(gcf, 'Color', 'w');
semilogy(f_gauge(valid_idx), psd_gauge(valid_idx), 'b-', 'LineWidth', 1.8);

xlabel('Frequency [Hz]');
ylabel('PSD [m^2/Hz]');
title('Power Spectral Density (Log Y) - Gauge 18 (Set 1)');
grid on;
grid minor;
xlim([0 freq_limit]);
ylim([1e-8 max(psd_gauge(valid_idx))*1.2]);
legend('Gauge 18');

%%
%% Amplitude Spectrum (Eta-f) for Gauge 18 (Set 1)

% 数据
eta_signal = eta_gauges_ds{1}(:,18); % Set 1 中 Gauge 18
N = length(eta_signal);
fs = SR; % Radar采样率（10Hz）
t = (0:N-1) / fs;

% 快速傅里叶变换（FFT）
Y = fft(eta_signal);
f = fs * (0:floor(N/2)) / N; % 单边频率轴

% 单边幅度谱
amp_spectrum = abs(Y(1:floor(N/2)+1)) / N;
amp_spectrum(2:end-1) = 2 * amp_spectrum(2:end-1);

% 限制频率到0-4.5 Hz
freq_limit = 2;
valid_idx = f <= freq_limit;

% 绘制
figure;
set(gcf, 'Color', 'w');
plot(f(valid_idx), amp_spectrum(valid_idx), 'b-', 'LineWidth', 1.8);

xlabel('Frequency [Hz]');
ylabel('Amplitude [m]');
title('Amplitude Spectrum (Eta-f) - Gauge 18 (Set 1)');
grid on;
grid minor;
xlim([0 freq_limit]);
ylim([0 max(amp_spectrum(valid_idx))*1.1]);
legend('Gauge 18');

%%
%% Amplitude Spectrum (Eta-f) - Log Y Axis

% 数据
eta_signal = eta_gauges_ds{1}(:,18); % Set 1 中 Gauge 18
N = length(eta_signal);
fs = SR; % Radar采样率（10Hz）
t = (0:N-1) / fs;

% FFT
Y = fft(eta_signal);
f = fs * (0:floor(N/2)) / N; % 单边频率轴

% 单边幅度谱
amp_spectrum = abs(Y(1:floor(N/2)+1)) / N;
amp_spectrum(2:end-1) = 2 * amp_spectrum(2:end-1);

% 限制频率到0-4.5 Hz
freq_limit = 4.5;
valid_idx = f <= freq_limit;

% 绘图
figure;
set(gcf, 'Color', 'w');
semilogy(f(valid_idx), amp_spectrum(valid_idx), 'b-', 'LineWidth', 1.8);

xlabel('Frequency [Hz]');
ylabel('Amplitude [m]');
title('Amplitude Spectrum (Eta-f) - Gauge 18 (Set 1) (Log Y)');
grid on;
grid minor;
xlim([0 freq_limit]);
ylim([1e-5, max(amp_spectrum(valid_idx))*1.2]);
legend('Gauge 18');
%%