1 %{
  2 mergePxxPlot.m
  3 
  4 This script processes and visualizes pre-multiplied power spectral density (PSD) data 
  5 from a merged dataset. The script performs the following tasks:
  6 
  7 1. **Load Data**:
  8     - Loads the merged PSD data from a `.mat` file (`merged_pxx.mat`).
  9 
 10 2. **Preprocess Data**:
 11     - Extracts frequency (`f_merge`), vertical positions (`Y_merge`), and PSD values (`pxx_merge`).
 12     - Computes pre-multiplied PSD values (`f * S_uu(f)`).
 13 
 14 3. **Generate Individual PSD Plots**:
 15     - Creates and saves individual PSD plots for each vertical position (`z`).
 16     - Saves the plots in a specified output folder (`mergePSD-figure`) as `.png` files.
 17 
 18 4. **Generate Combined PSD Plot**:
 19     - Selects specific vertical positions (`z`) for visualization.
 20     - Normalizes the selected `z` values for color mapping.
 21     - Plots the selected PSD data with a color gradient based on vertical position.
 22     - Adds a colorbar to indicate the mapping of colors to `z` values.
 23     - Saves the combined plot as `.fig`, `.emf`, and high-resolution `.jpg` files.
 24 
 25 5. **Interpolate Data**:
 26     - Interpolates the PSD data onto a finer grid for contour plotting.
 27     - Uses logarithmic spacing for frequency (`f`) and linear spacing for vertical positions (`z`).
 28 
 29 6. **Generate Contour Plot**:
 30     - Creates a filled contour plot of the interpolated PSD data.
 31     - Configures axes, labels, and colorbar with LaTeX formatting.
 32     - Uses a custom colormap (`sky`) for visualization.
 33 
 34 Dependencies:
 35 - Requires the `merged_pxx.mat` file containing the following variables:
 36   - `f_merge`: Frequency data (Hz).
 37   - `Y_merge`: Vertical positions (m).
 38   - `pxx_merge`: PSD values.
 39 - Requires the `sky` colormap function for custom color mapping.
 40 
 41 Output:
 42 - Individual PSD plots saved in the `mergePSD-figure` folder.
 43 - Combined PSD plot saved as `.fig`, `.emf`, and `.jpg` files.
 44 - Contour plot displayed with interpolated PSD data.
 45 
 46 Notes:
 47 - Ensure the `merged_pxx.mat` file is in the same directory as the script or provide the correct path.
 48 - The `sky` colormap function must be available in the MATLAB path.
 49 
 50 Author: L-N1988
 51 Date: 2025-3-30
 52 %}
 53 clc; clear; close all;
 54 mergeData = load('merged_pxx.mat');
 55 
 56 % Measuring points along vertical centre line
 57 xv = mergeData.f_merge;
 58 yv = repmat(double(mergeData.Y_merge), 1, size(xv, 2));
 59 vv = xv .* mergeData.pxx_merge; % pre-PSD
 60 
 61 % Plot and save the merged PSD
 62 outputFolder = 'mergePSD-figure';
 63 if ~exist(outputFolder, 'dir')
 64     mkdir(outputFolder);
 65 end
 66 
 67 for ii = 1:length(mergeData.Y_merge)
 68     figure;
 69     plot(xv(ii, :), vv(ii, :), 'LineWidth', 2);
 70     set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
 71     set(gca, 'FontSize', 16);
 72     set(xlabel("$f$ (Hz)"), 'Interpreter', 'latex');
 73     set(ylabel("$fS_{uu}(f) (\rm m^2/s^2)$"), 'Interpreter', 'latex');
 74     set(title(sprintf("PSD at $z = %.3f$ m", mergeData.Y_merge(ii))), 'Interpreter', 'latex');
 75     saveas(gcf, fullfile(outputFolder, sprintf("merged_psd_z%.3f.png", mergeData.Y_merge(ii))));
 76 end
 77 
 78 %%
 79 % Create a single figure
 80 figure;
 81 
 82 % Define the indices for the five data series
 83 indices = 10:10:length(mergeData.Y_merge); % Select every 5th index
 84 
 85 % Extract the corresponding z values for these indices
 86 z_values = mergeData.Y_merge(indices); % Vertical positions for the selected indices
 87 
 88 % Normalize the selected z values to [0, 1] for color mapping
 89 z_min = min(z_values);
 90 z_max = max(z_values);
 91 z_normalized = (z_values - z_min) / (z_max - z_min); % Normalize to [0, 1]
 92 
 93 % Define a colormap (e.g., 'jet') for the five selected series
 94 cmap = sky(length(indices)); % Colormap with 5 colors (one for each selected series)
 95 
 96 % Plot the selected data series with a color gradient
 97 hold on;
 98 for idx = 1:length(indices)
 99     ii = indices(idx); % Get the actual index
100     % Get the color for this vertical position
101     color_idx = round(z_normalized(idx) * (size(cmap, 1) - 1)) + 1; % Map to colormap index
102     plot(xv(ii, :), vv(ii, :), 'LineWidth', 1.5, 'Color', cmap(color_idx, :));
103 end
104 hold off;
105 
106 % Set the axes properties
107 set(gca, 'XScale', 'log');
108 set(gca, 'YScale', 'log');
109 set(gca, 'FontSize', 16);
110 
111 % Set labels and title with LaTeX interpreter
112 set(xlabel('$f$ (Hz)'), 'Interpreter', 'latex', 'FontSize', 16);
113 set(ylabel('$fS_{uu}(f) (\rm m^2/s^2)$'), 'Interpreter', 'latex', 'FontSize', 16);
114 set(title('Pre-multiplied Power Spectral Density'), 'Interpreter', 'latex', 'FontSize', 16);
115 
116 % Add a colorbar to show the mapping of colors to z values
117 colormap(cmap);
118 cbar = colorbar;
119 set(cbar, 'FontSize', 16);
120 set(get(cbar, 'Label'), 'String', '$z$ (m)', 'Interpreter', 'latex', 'FontSize', 16);
121 clim([z_min z_max]); % Set the colorbar limits to the range of selected z values
122 
123 % Save the figure
124 base_filename = fullfile(outputFolder, 'merged_psd_selected');
125 % Save as MATLAB figure (.fig)
126 savefig(gcf, [base_filename '.fig']);
127 
128 % Save as EMF (.emf)
129 print(gcf, [base_filename '.emf'], '-dmeta');
130 
131 % Save as high-DPI JPEG using print (e.g., 500 DPI)
132 print(gcf, [base_filename '.jpg'], '-djpeg', '-r500');
133 %%
134 
135 % Interpolate data
136 [grid_row, grid_col] = deal(max(400, size(yv, 1)*10), max(50, round(size(xv, 2)/100)));
137 xq = logspace(...
138     log10(min(xv(xv > 0))), log10(max(xv(:))), ...
139     grid_col);
140 yq = linspace(...
141     min(yv(:)), max(yv(:)), ...
142     grid_row);
143 [xq, yq] = meshgrid(xq, yq);
144 
145 vq = griddata(xv, yv, vv, xq, yq);
146 
147 %% Plot contour
148 contourf(xq, yq, vq, 10, 'LineStyle', '--');
149 set(gca, 'XScale', 'log'); set(gca, 'FontSize', 16); %set(gca, 'YScale', 'log');
150 set(xlabel("$f$ (Hz)"), 'Interpreter', 'latex');
151 set(ylabel("$z(\rm m)$"), 'Interpreter', 'latex');
152 colormap("sky");
153 col = colorbar();
154 set(ylabel(col,"$fS_{uu}(f) (\rm m^2/s^2)$"), 'Interpreter', 'latex');