A Matlab code of how to do dynamic mode decomposition of a video sequence. (ChatGPT)
% DMD on a Video Sequence clear; close all; clc;
%% Load Video videoFile = 'video.mp4'; % Change to your video file v = VideoReader(videoFile);
% Parameters resizeFactor = 0.25; % Reduce size to speed up numFrames = min(100, floor(v.Duration * v.FrameRate)); % Limit number of frames
% Read and preprocess video frames X = []; for i = 1:numFrames frame = readFrame(v); grayFrame = im2double(rgb2gray(frame)); % Convert to grayscale grayFrame = imresize(grayFrame, resizeFactor); % Resize X(:,i) = grayFrame(:); % Vectorize and store as column end
%% Build data matrices X1 and X2 X1 = X(:, 1:end-1); X2 = X(:, 2:end);
%% Singular Value Decomposition [U, S, V] = svd(X1, 'econ');
% Choose r - number of modes r = 20; % Adjust based on singular value decay U_r = U(:, 1:r); S_r = S(1:r, 1:r); V_r = V(:, 1:r);
P1. Can you give me a Matlab code for Dynamic Mode Decomposition ? The code should be nicely commented so one can understand the purpose of each section.
P2. Modify the Matlab code to handle an array of flattened image data
P3. For the video sequence, how many snapshots, captured at what fps will be needed to show a mode at 200 Hz and a mode at 2000 Hz?
P4. So DMD gives the eigen modes, their amplitudes and their complex frequencies. From the complex frequencies, the real part gives the growth or decay rate and the imaginary part gives the circular frequncy. Is this interpretation correct ?
P5. How to get the amplitude of a mode ?
P6. Can you give an example of complex amplitude ?
The Bychkov limit, also known as the detonation cell size limit, is a concept in detonation physics that describes the minimum size of a detonation cell needed to sustain a propagating detonation wave. This limit is not fixed but depends on factors like the mixture composition, pressure, and tube diameter. [1, 1, 2, 2, 3, 3, 4, 4]. Detonation Cell Size: In a detonation, the reaction zone behind the shock wave propagates in a cellular structure, with each cell representing a single, self-sustaining detonation event. [1, 1, 3, 3, 5, 6]. Bychkov Limit: The Bychkov limit is a critical length scale related to the detonation cell size. It represents the minimum size of a detonation cell that can propagate a stable and sustained detonation wave. [1, 1, 2, 2, 7, 8]. Factors Affecting the Limit: [1, 2, 3, 9] (a) Mixture Composition: The type and ratio of fuel and oxidizer significantly impact the detonation cell size and thus the Bychkov limit. [1, 1, 3, 3, 9, 9] (b) Pressure: Higher pressures generally lead to smaller detonation cell sizes and a lower Bychkov limit. [1, 1, 3, 3, 9, 9, 10, 11, 12] (c) Tube Diameter: The diameter of the container or tube where the detonation is propagating can influence the minimum cell size required for sustained detonation. [1, 1, 2, 2, 3, 3]. Importance: Understanding the Bychkov limit is crucial for designing and analyzing detonation processes, predicting detonability, and mitigating the risks associated with explosions. [2, 2, 3, 3, 4, 4] (a) Experimental and Numerical Studies: The Bychkov limit and its relation to detonation cell size are extensively studied through experiments (e.g., in detonation tubes) and numerical simulations. [13, 13, 14, 14, 15, 15] (b) Examples: In experiments with hydrogen-oxygen mixtures, the Bychkov limit has been observed to be in the range of a few millimeters to several centimeters, depending on the specific conditions. [13, 13, 15, 15, 16, 16] (c) Transition to Deflagration: When the detonation cell size exceeds the Bychkov limit, the detonation wave may become unstable and transition to a deflagration (a slower, subsonic combustion process). [13, 13, 17, 17, 18, 18, 19]
%Here's a MATLAB code that performs Independent Component Analysis (ICA) on a video sequence. The goal is to treat each frame as a data %vector and then extract statistically independent components from this set using ICA.
%This is particularly useful for motion segmentation, source separation, or background subtraction tasks.
% ICA on a video sequence clc; clear; close all;
% Load video videoFile = 'your_video.mp4'; % Change to your video file v = VideoReader(videoFile); numFrames = floor(v.Duration * v.FrameRate);
% Read first frame to get dimensions frame = readFrame(v); grayFrame = rgb2gray(frame); [rows, cols] = size(grayFrame); frameSize = rows * cols;
% Reset video to the beginning v.CurrentTime = 0;
% Construct data matrix: each column is a vectorized frame dataMatrix = zeros(frameSize, numFrames); frameCount = 0;
while hasFrame(v) frame = readFrame(v); grayFrame = rgb2gray(frame); grayFrame = im2double(grayFrame); % Normalize to [0,1] frameCount = frameCount 1; dataMatrix(:, frameCount) = grayFrame(:); end
% Remove mean for ICA X = dataMatrix - mean(dataMatrix, 2);
% Perform ICA using FastICA addpath('fastica'); % Make sure FastICA is in your path (download from https://research.ics.aalto.fi/ica/fastica/) [icasig, A, W] = fastica(X, 'numOfIC', 10); % You can change the number of ICs
% Display some ICs as images figure; for i = 1:9 subplot(3,3,i); imagesc(reshape(icasig(i, :), [1, frameCount])); colormap gray; axis off; title(['IC ', num2str(i)]); end
% You can also reshape the spatial ICs if desired figure; for i = 1:9 subplot(3,3,i); imagesc(reshape(A(:, i), rows, cols)); colormap gray; axis off; title(['Spatial IC ', num2str(i)]); end
% Set image directory imageDir = 'images/'; % Change this to your image folder path imageFiles = dir(fullfile(imageDir, '*.png')); % Change extension if needed (e.g., .jpg, .bmp) numImages = length(imageFiles);
% Read first image to get dimensions firstImage = imread(fullfile(imageDir, imageFiles(1).name)); if size(firstImage, 3) == 3 firstImage = rgb2gray(firstImage); % Convert to grayscale if RGB end [rows, cols] = size(firstImage); imageSize = rows * cols;
% Create data matrix where each column is an image vector dataMatrix = zeros(imageSize, numImages);
for i = 1:numImages img = imread(fullfile(imageDir, imageFiles(i).name)); if size(img, 3) == 3 img = rgb2gray(img); end img = im2double(img); % Convert to double dataMatrix(:, i) = img(:); % Vectorize image end
% Mean normalization meanImage = mean(dataMatrix, 2); X = dataMatrix - meanImage;
% Compute covariance matrix C = cov(X');
% Eigen decomposition [eigVectors, eigValues] = eig(C); eigValues = diag(eigValues);
% Sort eigenvalues and vectors in descending order [~, idx] = sort(eigValues, 'descend'); eigVectors = eigVectors(:, idx);
% Project data onto principal components k = 10; % Number of principal components to keep eigVectors_k = eigVectors(:, 1:k); projectedData = eigVectors_k' * X;
% Reconstruct an image using top k components reconstructed = eigVectors_k * projectedData meanImage;
% Show original and reconstructed images figure; subplot(1, 2, 1); imshow(reshape(dataMatrix(:, 1), rows, cols), []); title('Original Image');
Above Image is from Arts, L.P.A., van den Broek, E.L. The fast continuous wavelet transformation (fCWT) for real-time, high-quality, noise-resistant time–frequency analysis. Nat Comput Sci 2, 47–58 (2022). https://doi.org/10.1038/s43588-021-00183-z
How does wavelet transform differ from spectrogram when used on the same time series data ?
ChatGPT : Both wavelet transform and spectrograms are time-frequency analysis tools, but they differ significantly in how they analyze and represent the frequency content of time series data. Here's a breakdown of the differences:
1. Spectrogram (based on Short-Time Fourier Transform - STFT)
How it works: - Splits the signal into short, fixed-length time windows. - Applies the Fourier Transform to each window. - Gives a time vs frequency plot showing how frequency content changes over time.
Key Characteristics: - Fixed time-frequency resolution: The window length is constant, so you have a trade-off: - Good time resolution → poor frequency resolution. - Good frequency resolution → poor time resolution. - Best for stationary or slowly varying signals.
Analogy: Like taking snapshots of a signal using the same size lens—details get lost either in time or frequency depending on the window size.
How it works: - Decomposes the signal using scaled and shifted versions of a wavelet (a small waveform). - Automatically adjusts the resolution: - High frequencies → short time windows (good time resolution, poor frequency resolution). - Low frequencies → long time windows (good frequency resolution, poor time resolution).
Key Characteristics: - Multi-resolution analysis: Adapts to different scales. - Better suited for non-stationary or transient signals. - Captures short-lived features like bursts, edges, or spikes better than STFT.
Analogy: Like zooming in and out on a signal with a flexible lens—details appear at multiple scales.
Comparison Table
| Feature | Spectrogram (STFT) | Wavelet Transform (CWT) | |--------|---------------------|-------------------------| | Resolution | Fixed | Adaptive (multi-resolution) | | Time-Frequency Tradeoff | Yes (fixed window) | No (scale-dependent) | | Frequency Analysis | Global (within each window) | Localized (scale-wise) | | Best for | Stationary signals | Non-stationary/transient signals | | Window Function | Fixed (e.g., Hamming, Hann) | Wavelet (e.g., Morlet, Mexican Hat) | | Output | Time-Frequency plot | Time-Scale plot (can be converted to frequency) |
In Practice
- Spectrograms are simpler and computationally cheaper, good for general overview. - Wavelet transforms give richer detail and are often used in applications like: - Biomedical signal analysis (e.g., ECG, EEG) - Seismic analysis - Fault detection in mechanical systems
If you're working with signals where the frequency content changes quickly (like engine vibrations, fault signatures, or speech), wavelets often reveal more subtle features than spectrograms.
Aerodynamics, Thermodynamics, Heat Transfer and Combustion Chemistry. The building blocks of propusion systems. While feedback control is also an integral part, courses like these primarily cover the mechanical aspects of propulsion systems.
1157 0 0 0
A Matlab code of how to do dynamic mode decomposition of a video sequence.
Amardip Ghosh #Advanced Propulsion Systems (APSYS) Lab
A Matlab code of how to do dynamic mode decomposition of a video sequence. (ChatGPT)
% DMD on a Video Sequence
clear; close all; clc;
%% Load Video
videoFile = 'video.mp4'; % Change to your video file
v = VideoReader(videoFile);
% Parameters
resizeFactor = 0.25; % Reduce size to speed up
numFrames = min(100, floor(v.Duration * v.FrameRate)); % Limit number of frames
% Read and preprocess video frames
X = [];
for i = 1:numFrames
frame = readFrame(v);
grayFrame = im2double(rgb2gray(frame)); % Convert to grayscale
grayFrame = imresize(grayFrame, resizeFactor); % Resize
X(:,i) = grayFrame(:); % Vectorize and store as column
end
%% Build data matrices X1 and X2
X1 = X(:, 1:end-1);
X2 = X(:, 2:end);
%% Singular Value Decomposition
[U, S, V] = svd(X1, 'econ');
% Choose r - number of modes
r = 20; % Adjust based on singular value decay
U_r = U(:, 1:r);
S_r = S(1:r, 1:r);
V_r = V(:, 1:r);
%% Build A tilde
A_tilde = U_r' * X2 * V_r / S_r;
%% Eigen decomposition
[W, D] = eig(A_tilde);
Phi = X2 * V_r / S_r * W; % DMD modes
%% Compute DMD eigenvalues and frequencies
lambda = diag(D);
omega = log(lambda);
%% Compute initial amplitudes b
x1 = X(:,1);
b = Phi x1;
%% Time dynamics
time_dynamics = zeros(r, numFrames - 1);
t = (0:numFrames-2) / v.FrameRate; % Time vector in seconds
for i = 1:length(t)
time_dynamics(:, i) = (b .* exp(omega * t(i)));
end
X_dmd = real(Phi * time_dynamics);
%% Visualize original vs reconstructed frame
originalFrame = reshape(X(:, 30), size(grayFrame));
reconstructedFrame = reshape(X_dmd(:, 30), size(grayFrame));
figure;
subplot(1,2,1); imshow(originalFrame,[]); title('Original Frame 30');
subplot(1,2,2); imshow(reconstructedFrame,[]); title('DMD Reconstructed Frame 30');
Some prompts I used for ChatGPT
P1. Can you give me a Matlab code for Dynamic Mode Decomposition ? The code should be nicely commented so one can understand the purpose of each section.
P2. Modify the Matlab code to handle an array of flattened image data
P3. For the video sequence, how many snapshots, captured at what fps will be needed to show a mode at 200 Hz and a mode at 2000 Hz?
P4. So DMD gives the eigen modes, their amplitudes and their complex frequencies. From the complex frequencies, the real part gives the growth or decay rate and the imaginary part gives the circular frequncy. Is this interpretation correct ?
P5. How to get the amplitude of a mode ?
P6. Can you give an example of complex amplitude ?
888 0 0 0
Bychkov Limit (ChatGPT, Google AI)
Amardip Ghosh #Advanced Propulsion Systems (APSYS) Lab
977 0 0 0
Independent Component Analysis (ICA) for a video sequence. ChatGPT
Amardip Ghosh #Advanced Propulsion Systems (APSYS) Lab
%Here's a MATLAB code that performs Independent Component Analysis (ICA) on a video sequence. The goal is to treat each frame as a data %vector and then extract statistically independent components from this set using ICA.
%This is particularly useful for motion segmentation, source separation, or background subtraction tasks.
% ICA on a video sequence
clc; clear; close all;
% Load video
videoFile = 'your_video.mp4'; % Change to your video file
v = VideoReader(videoFile);
numFrames = floor(v.Duration * v.FrameRate);
% Read first frame to get dimensions
frame = readFrame(v);
grayFrame = rgb2gray(frame);
[rows, cols] = size(grayFrame);
frameSize = rows * cols;
% Reset video to the beginning
v.CurrentTime = 0;
% Construct data matrix: each column is a vectorized frame
dataMatrix = zeros(frameSize, numFrames);
frameCount = 0;
while hasFrame(v)
frame = readFrame(v);
grayFrame = rgb2gray(frame);
grayFrame = im2double(grayFrame); % Normalize to [0,1]
frameCount = frameCount 1;
dataMatrix(:, frameCount) = grayFrame(:);
end
% Remove mean for ICA
X = dataMatrix - mean(dataMatrix, 2);
% Perform ICA using FastICA
addpath('fastica'); % Make sure FastICA is in your path (download from https://research.ics.aalto.fi/ica/fastica/)
[icasig, A, W] = fastica(X, 'numOfIC', 10); % You can change the number of ICs
% Display some ICs as images
figure;
for i = 1:9
subplot(3,3,i);
imagesc(reshape(icasig(i, :), [1, frameCount]));
colormap gray; axis off;
title(['IC ', num2str(i)]);
end
% You can also reshape the spatial ICs if desired
figure;
for i = 1:9
subplot(3,3,i);
imagesc(reshape(A(:, i), rows, cols));
colormap gray; axis off;
title(['Spatial IC ', num2str(i)]);
end
disp('ICA completed successfully!');
992 0 0 0
Principle Component Analysis (PCA) Matlab Code (ChatGPT)
Amardip Ghosh #Advanced Propulsion Systems (APSYS) Lab
% PCA on an image set
clc; clear; close all;
% Set image directory
imageDir = 'images/'; % Change this to your image folder path
imageFiles = dir(fullfile(imageDir, '*.png')); % Change extension if needed (e.g., .jpg, .bmp)
numImages = length(imageFiles);
% Read first image to get dimensions
firstImage = imread(fullfile(imageDir, imageFiles(1).name));
if size(firstImage, 3) == 3
firstImage = rgb2gray(firstImage); % Convert to grayscale if RGB
end
[rows, cols] = size(firstImage);
imageSize = rows * cols;
% Create data matrix where each column is an image vector
dataMatrix = zeros(imageSize, numImages);
for i = 1:numImages
img = imread(fullfile(imageDir, imageFiles(i).name));
if size(img, 3) == 3
img = rgb2gray(img);
end
img = im2double(img); % Convert to double
dataMatrix(:, i) = img(:); % Vectorize image
end
% Mean normalization
meanImage = mean(dataMatrix, 2);
X = dataMatrix - meanImage;
% Compute covariance matrix
C = cov(X');
% Eigen decomposition
[eigVectors, eigValues] = eig(C);
eigValues = diag(eigValues);
% Sort eigenvalues and vectors in descending order
[~, idx] = sort(eigValues, 'descend');
eigVectors = eigVectors(:, idx);
% Project data onto principal components
k = 10; % Number of principal components to keep
eigVectors_k = eigVectors(:, 1:k);
projectedData = eigVectors_k' * X;
% Reconstruct an image using top k components
reconstructed = eigVectors_k * projectedData meanImage;
% Show original and reconstructed images
figure;
subplot(1, 2, 1);
imshow(reshape(dataMatrix(:, 1), rows, cols), []);
title('Original Image');
subplot(1, 2, 2);
imshow(reshape(reconstructed(:, 1), rows, cols), []);
title(['Reconstructed with ', num2str(k), ' PCs']);
disp('PCA completed successfully!');
1137 0 0 0
Wavelets in Signal Processing
Amardip Ghosh #Advanced Propulsion Systems (APSYS) Lab
Above Image is from Arts, L.P.A., van den Broek, E.L. The fast continuous wavelet transformation (fCWT) for real-time, high-quality, noise-resistant time–frequency analysis. Nat Comput Sci 2, 47–58 (2022). https://doi.org/10.1038/s43588-021-00183-z
How does wavelet transform differ from spectrogram when used on the same time series data ?
ChatGPT : Both wavelet transform and spectrograms are time-frequency analysis tools, but they differ significantly in how they analyze and represent the frequency content of time series data. Here's a breakdown of the differences:
1. Spectrogram (based on Short-Time Fourier Transform - STFT)
How it works:
- Splits the signal into short, fixed-length time windows.
- Applies the Fourier Transform to each window.
- Gives a time vs frequency plot showing how frequency content changes over time.
Key Characteristics:
- Fixed time-frequency resolution: The window length is constant, so you have a trade-off:
- Good time resolution → poor frequency resolution.
- Good frequency resolution → poor time resolution.
- Best for stationary or slowly varying signals.
Analogy: Like taking snapshots of a signal using the same size lens—details get lost either in time or frequency depending on the window size.
2. Wavelet Transform (especially Continuous Wavelet Transform - CWT)
How it works:
- Decomposes the signal using scaled and shifted versions of a wavelet (a small waveform).
- Automatically adjusts the resolution:
- High frequencies → short time windows (good time resolution, poor frequency resolution).
- Low frequencies → long time windows (good frequency resolution, poor time resolution).
Key Characteristics:
- Multi-resolution analysis: Adapts to different scales.
- Better suited for non-stationary or transient signals.
- Captures short-lived features like bursts, edges, or spikes better than STFT.
Analogy: Like zooming in and out on a signal with a flexible lens—details appear at multiple scales.
Comparison Table
| Feature | Spectrogram (STFT) | Wavelet Transform (CWT) |
|--------|---------------------|-------------------------|
| Resolution | Fixed | Adaptive (multi-resolution) |
| Time-Frequency Tradeoff | Yes (fixed window) | No (scale-dependent) |
| Frequency Analysis | Global (within each window) | Localized (scale-wise) |
| Best for | Stationary signals | Non-stationary/transient signals |
| Window Function | Fixed (e.g., Hamming, Hann) | Wavelet (e.g., Morlet, Mexican Hat) |
| Output | Time-Frequency plot | Time-Scale plot (can be converted to frequency) |
In Practice
- Spectrograms are simpler and computationally cheaper, good for general overview.
- Wavelet transforms give richer detail and are often used in applications like:
- Biomedical signal analysis (e.g., ECG, EEG)
- Seismic analysis
- Fault detection in mechanical systems
If you're working with signals where the frequency content changes quickly (like engine vibrations, fault signatures, or speech), wavelets often reveal more subtle features than spectrograms.
5834 4 3 3
Thermodynamics and Aerospace Propulsion Systems Class
Amardip Ghosh #Advanced Propulsion Systems (APSYS) Lab
Aerodynamics, Thermodynamics, Heat Transfer and Combustion Chemistry. The building blocks of propusion systems. While feedback control is also an integral part, courses like these primarily cover the mechanical aspects of propulsion systems.