1 % PSD.m
 2 %
 3 % This script processes velocity data to compute the Power Spectral Density (PSD)
 4 % using Welch's method. It performs the following steps:
 5 % 1. Clears the workspace, closes all figures, and clears the command window.
 6 % 2. Loads velocity data from a specified file.
 7 % 3. Computes the PSD for each spatial point in the velocity data using Welch's method.
 8 % 4. Applies a 3x3 moving average filter to smooth the PSD data.
 9 % 5. Saves the processed PSD data to a .mat file for further analysis.
10 %
11 % Key Variables:
12 % - casePath: Path to the case directory containing the input data.
13 % - Fs: Sampling frequency (Hz).
14 % - X, Y: Spatial coordinates of the velocity data.
15 % - xmesh, ymesh: Mesh grids for the spatial coordinates.
16 % - U_xt: Mean velocity over time.
17 % - u_pri: Primary velocity data (3D matrix).
18 % - window_length: Length of the window used in Welch's method.
19 % - nfft: Number of FFT points (next power of 2 of window_length).
20 % - fs: Frequency matrix (3D matrix).
21 % - pxxs: PSD matrix (3D matrix).
22 %
23 % Functions:
24 % - patch_moving_average_conv: Computes a moving average using convolution
25 %   with a specified kernel. Used to smooth the PSD data.
26 %
27 % Outputs:
28 % - Saves the processed PSD data (fs, pxxs) along with spatial coordinates
29 %   (X, Y, xmesh, ymesh) and mean velocity (U_xt) to 'pxxs.mat'.
30 %
31 % Notes:
32 % - The script uses a 3x3 moving average filter for smoothing.
33 % - Welch's method is applied to each spatial point in the velocity data.
34 % - The PSD is visualized on a log-log scale during computation.
35 % Clear workspace, close all figures, and clear command window
36 clc; clear; close all;
37 
38 % Define the path to the case directory and sampling frequency
39 casePath = '..';
40 Fs = 1100; % Sampling frequency (Hz) edited by YZ
41 
42 % Load the data from the specified file
43 data = load(fullfile(casePath, 'figure_data', 'u4pxx.mat'));
44 X = data.X; % X-coordinates
45 Y = data.Y; % Y-coordinates
46 xmesh = data.xmesh; % X-mesh grid
47 ymesh = data.ymesh; % Y-mesh grid
48 U_xt = mean(data.U_t, 2); % Mean velocity over time
49 u_pri = data.u_pri; % Primary velocity data
50 
51 % Get the dimensions of the velocity data
52 [m, n, ~] = size(u_pri);
53 
54 % Define the window length for Welch's method and calculate FFT parameters
55 window_length = 8000; % Length of the window for Welch's method
56 nfft = 2^nextpow2(window_length); % Next power of 2 for FFT
57 output_length = nfft/2 + 1; % Output length for real-valued input
58 
59 % Initialize matrices to store frequency and power spectral density (PSD)
60 fs = zeros(m, n, output_length); % Frequency matrix
61 pxxs = zeros(m, n, output_length); % PSD matrix
62 
63 % Loop through each spatial point to compute PSD using Welch's method
64 for ii = 1:m
65     for jj = 1:n
66         % Compute PSD and frequency using Welch's method
67         [pxxs(ii, jj, :), fs(ii, jj, :)] = pwelch(squeeze(u_pri(ii, jj, :)), ...
68             window_length, [], [], Fs);
69 
70         % Plot the PSD on a log-log scale for visualization
71         loglog(squeeze(fs(ii, jj, :)), squeeze(pxxs(ii, jj, :)));
72     end
73 end
74 
75 % Define the window size for moving average (3x3)
76 window_size = 3;
77 kernel = ones(window_size, window_size); % 3x3 matrix of ones
78 
79 % Apply a 3x3 moving average filter to smooth the PSD data
80 pxxs = patch_moving_average_conv(pxxs, kernel);
81 
82 % Save the processed data to a .mat file
83 save(fullfile(casePath, 'figure_data', 'pxxs.mat'), 'X', 'xmesh', 'Y', 'ymesh', 'fs', 'pxxs', 'U_xt');
84 
85 % Function to compute a moving average using convolution
86 function moving_avg_matrix = patch_moving_average_conv(matrix, kernel)
87     % Input: matrix (e.g., 6x10), kernel (e.g., 3x3)
88     % Output: moving_avg_matrix with the same size as input, containing averages
89 
90     % Perform convolution with 'same' to keep output size equal to input
91     conv_result = convn(matrix, kernel, 'same');
92 
93     % Create a matrix to count how many elements contribute at each position
94     contrib_count = convn(ones(size(matrix)), kernel, 'same');
95 
96     % Compute the moving average by dividing the sum by the number of contributing elements
97     moving_avg_matrix = conv_result ./ contrib_count;
98 end