1 %{
  2 mergePSD.m - A MATLAB script for merging and aligning power spectral density (PSD) data from two datasets.
  3 
  4 It demonstrates how to:
  5 1. Load and process PSD data from two sources.
  6 2. Align datasets based on their Y-axis values.
  7 3. Merge low-frequency and high-frequency data into a single dataset.
  8 
  9 Key Features:
 10 - Handles datasets with different lengths by finding the best alignment offset.
 11 - Slices and merges data based on a specified frequency cutoff.
 12 - Combines data while preserving the integrity of the frequency and PSD values.
 13 
 14 Functions:
 15 1. mergeAlighment:
 16     - Aligns and merges two datasets (low-frequency and high-frequency).
 17     - Slices data based on a frequency cutoff and aligns them in the Y direction.
 18     - Outputs merged frequency, PSD, and Y-axis data.
 19 
 20 Inputs:
 21 - Two datasets containing frequency (`f`), PSD (`pxx`), and Y-axis (`Y`) values.
 22 - A frequency cutoff (`cut_f`) to separate low and high-frequency data.
 23 
 24 Outputs:
 25 - Merged frequency (`f_merge`), PSD (`pxx_merge`), and Y-axis (`Y_merge`) data.
 26 
 27 Usage:
 28 - This script is useful for combining PSD data from different measurement ranges or resolutions.
 29 - It can be adapted for other types of data alignment and merging tasks.
 30 
 31 Note:
 32 - Ensure the input data files (`pxxs.mat`) contain the required variables (`pxxs`, `fs`, `Y`).
 33 - The script assumes the datasets are stored in specific directories (`C/figure_data` and `L/figure_data`).
 34 
 35 %}
 36 
 37 clc; clear; close all;
 38 data_c = load("C/figure_data/pxxs.mat");
 39 data_l = load("L/figure_data/pxxs.mat");
 40 
 41 pxx_c = data_c.pxxs; f_c = data_c.fs;
 42 % Find the middle index of the data
 43 mid_c = round(size(pxx_c, 2)/2);
 44 pxx_c = squeeze(pxx_c(:, mid_c, :));
 45 f_c = squeeze(f_c(:, mid_c, :));
 46 % Remove high frequency white noise data
 47 out_f = 100; % cut edge of frequency
 48 out_index = f_c(1, :) > out_f;
 49 f_c(:, out_index) = [];
 50 pxx_c(:, out_index) = [];
 51 
 52 pxx_l = data_l.pxxs; f_l = data_l.fs;
 53 mid_l = round(size(pxx_l, 2)/2);
 54 pxx_l = squeeze(pxx_l(:, mid_l, :));
 55 f_l = squeeze(f_l(:, mid_l, :));
 56 
 57 cut_f = 2; % cut edge of merged frequency
 58 index_c = (f_c >= cut_f);
 59 index_l = (f_l <= cut_f);
 60 
 61 Y_c = data_c.Y;
 62 Y_l = data_l.Y;
 63 
 64 % Align datasets based on Y-axis values
 65 if length(Y_l) > length(Y_c)
 66     longer = Y_l;
 67     shorter = Y_c;
 68 else
 69     longer = Y_c;
 70     shorter = Y_l;
 71 end
 72 best_offset = 0;
 73 min_diff = inf;
 74 for offset = 0:(length(longer) - length(shorter))
 75     current_diff = sum(abs(longer(1+offset:offset+length(shorter)) - shorter));
 76     if current_diff < min_diff
 77         best_offset = offset;
 78         min_diff = current_diff;
 79     end
 80 end
 81 
 82 if length(Y_l) > length(Y_c)
 83     [f_merge, pxx_merge, Y_merge] = mergeAlighment({f_l, pxx_l, Y_l}, {f_c, pxx_c, Y_c}, cut_f, best_offset);
 84 else
 85     [f_merge, pxx_merge, Y_merge] = mergeAlighment({f_c, pxx_c, Y_c}, {f_l, pxx_l, Y_l}, cut_f, best_offset);
 86 end
 87 
 88 save('merged_pxx.mat', 'f_merge', 'pxx_merge', 'Y_merge');
 89 
 90 function [f_merge, pxx_merge, Y_merge] = mergeAlighment(longer, shorter, cut_f, offset)
 91     if max(longer{1}, [], "all") < max(shorter{1}, [], "all")
 92         low_freq_data = longer;
 93         high_freq_data = shorter;
 94         [f_l, pxx_l, Y_l] = low_freq_data{:};
 95         [f_c, pxx_c, Y_c] = high_freq_data{:};
 96         % Slice low frequency data to align with high frequency data in Y direction
 97         pxx_l = pxx_l(offset+1:offset+length(Y_c), :);
 98         f_l = f_l(offset+1:offset+length(Y_c), :);
 99         Y_l = Y_l(offset+1:offset+length(Y_c), :);
100     else
101         low_freq_data = shorter;
102         high_freq_data = longer;
103         [f_l, pxx_l, Y_l] = low_freq_data{:};
104         [f_c, pxx_c, Y_c] = high_freq_data{:};
105         % Slice high frequency data to align with low frequency data in Y direction
106         pxx_c = pxx_c(offset+1:offset+length(Y_l), :);
107         f_c = f_c(offset+1:offset+length(Y_l), :);
108         Y_c = Y_c(offset+1:offset+length(Y_l), :);
109     end
110     Y_merge = mean([Y_l, Y_c], 2);
111 
112     % Preallocate memory for f_merge and pxx_merge
113     max_columns = size(f_l, 2) + size(f_c, 2); % Maximum possible columns
114     f_merge = zeros(length(Y_merge), max_columns);
115     pxx_merge = zeros(length(Y_merge), max_columns);
116 
117     % Slice low frequency data based on cut_f
118     for ii = 1:length(Y_merge)
119         index_l = (f_l(ii, :) <= cut_f);
120         index_c = (f_c(ii, :) >= cut_f);
121         merged_f = [f_l(ii, index_l) f_c(ii, index_c)];
122         merged_pxx = [pxx_l(ii, index_l) pxx_c(ii, index_c)];
123         assert(length(merged_f) == length(merged_pxx), "Length mismatch between frequency and PSD data.");
124 
125         % Store merged data
126         f_merge(ii, 1:length(merged_f)) = merged_f;
127         pxx_merge(ii, 1:length(merged_pxx)) = merged_pxx;
128         plot(f_c(ii, :), pxx_c(ii, :));
129         hold on;
130         plot(f_l(ii, :), pxx_l(ii, :));
131         plot(merged_f, merged_pxx);
132         set(gca, 'XScale', 'log', 'YScale', 'log'); hold off;
133     end
134 
135     % Trim unused columns
136     f_merge = f_merge(:, 1:length(merged_f));
137     pxx_merge = pxx_merge(:, 1:length(merged_pxx));
138 end