Do you have access to MATLAB or Octave or similar?
In general, I would:
- Compute the spectrogram (STFT) of the signal you're interested in.
- Perform any required scaling to get a good amplitude estimate.
- Apply some thresholding to treat the noise floor.
- Devise a scheme for identifying the amplitude maxima of the harmonics you're interested in.
This final point will be much easier if you know the fundamental. For example, if you're playing a note at 440Hz (A440), and you know you only want to examine the fundamental, second, and third harmonics, you can easily define some brackets to search within.
If you don't know the fundamental, you will have to implement some method for determining what the fundamental is. However, assuming you do know the fundamental, I've provided a MATLAB example.
This example first builds an arbitrary signal with some sinusoids that decay at differing rates, then analyses the signal to identify harmonic amplitudes over time.
%% Generate a signal with harmonics.
% Build fundamental and harmonic signals.
samplePeriod = 0.001;
time = 0:samplePeriod:5; % [s]
harmonics = [1, 3, 5, 7];
fundamentalFreq = 50; % [Hz]
componentFreqs = fundamentalFreq .* harmonics; % [Hz]
componentFreqs = componentFreqs .* 2 .* pi; % [rad/s]
componentFreqs = repmat(componentFreqs', 1, numel(time));
componentSignals = sin(componentFreqs .* repmat(time, numel(harmonics), 1));
% Apply decay to the harmonics.
decay = [0, 0.1, 0.5, 1];
harmonicAmplitude = zeros(numel(harmonics), numel(time));
for iHarmonic = 1 : numel(harmonics)
harmonicAmplitude(iHarmonic, :) = exp(-1 * decay(iHarmonic) * time);
componentSignals(iHarmonic, :) = componentSignals(iHarmonic, :) .* harmonicAmplitude(iHarmonic, :) .* 1;
end
% Generate final signal.
signal = sum(componentSignals, 1);
%% Perform time-frequency analysis.
sampleFrequency = 1 / samplePeriod;
% Compute the spectrogram.
NFFT = 2^12;
winDuration = 0.1; % [s]
winLength = winDuration * sampleFrequency;
winFactor = mean(hamming(winLength));
[S, freq, specTime, p] = spectrogram(signal, winLength, [], NFFT, sampleFrequency);
% Compute FFT amplitude.
amplitude = abs(S) .* 2 ./ winLength ./ winFactor;
% Threshold.
floorLevel = max(max(amplitude)) * 0.1;
amplitude(amplitude < floorLevel) = 0;
% Determine the amplitude of the harmonics of interest.
fundamentalFreq = 50; % [Repeat of above]
harmonicsOfInterest = 1 : 10;
harmonicSearchWidth = 5; % [Hz]
harmonicBracketLower = harmonicsOfInterest .* fundamentalFreq - harmonicSearchWidth / 2;
harmonicBracketUpper = harmonicsOfInterest .* fundamentalFreq + harmonicSearchWidth / 2;
% Find the amplitude maxima for each harmonic at each point in time.
estimatedHarmonicAmplitude = zeros(numel(specTime), numel(harmonicsOfInterest));
for iHarmonic = 1 : numel(harmonicsOfInterest)
activeBracket = freq >= harmonicBracketLower(iHarmonic) & freq <= harmonicBracketUpper(iHarmonic);
bracketAmplitude = amplitude(activeBracket, :);
estimatedHarmonicAmplitude(:, iHarmonic) = max(bracketAmplitude, [], 1)';
end
% Visualise output.
figure;
surf(specTime, freq, amplitude, 'EdgeColor', 'none');
xlabel('Time [s]');
ylabel('Frequency [Hz]');
zlabel('Amplitude');
view(2);
figure;
plot(specTime, estimatedHarmonicAmplitude);
xlabel('Time [s]');
ylabel('Amplitude');
box on;
grid on;
legendBuilder = @(a) ['Harmonic: ', num2str(a, '%1.0f')];
for iHarmonic = 1 : numel(harmonicsOfInterest)
legendString{iHarmonic} = legendBuilder(harmonicsOfInterest(iHarmonic));
end
legend(legendString);
