I'll propose a short answer and develop it if you need more information
So, step 1 and 3 are somehow easy to implement
I would map step 2 into:
2.1) From each wave extract a set of features. This step basically defines what it means for tow waves to be similar and what information from waves is irrelevant.
2.2) Cluster the waves in buckets. Choose a method that computes the number of clusters automatically (you still have to specify a value that roughly defines the similarity between the data)
For 2.1 try to use more than one algorithm that can describe your wave. There are a wide variety of algorithms and your FFT ideea is good (that is, for each wave, compute FFT and use first k coefficients as features). Note that some algorithms generate features that are phase invariant, other keep the phase information. The choice is yours, whenever the phase information is useful or not.
For 2.2 you can use out of the box clustering algorithms. I suggest to start with the simpler ones, than climb up to more sophisticated approaches. There are mainly two categories, one for which you state the desired number of clusters and the other one where you specify something equivalent of a threshold and the algorithm clusters the data in as many buckets as necessary. Start with a k-means clustering and stop maybe at a mean shift clustering. Hierarchical clustering is another idea but could be hard to work with it from the beginning.
From my perspective step 2.2 is the most difficult to get it right. The meta parameter that controls the clustering is tricky to be tuned. Ideally you should have a way to measure independently how good the clustering is (ex. make a hand made bucketing that will be a ground truth for your method or derive some domain specific measures for the images). Based on the ground truth you can "tune" the meta parameter to an optimal value. The tuning is done using cross-validation or a simple train/evaluate splitting.
p.s. The clustering algorithm basically computes some distances between the features. You can measure directly the distance between two 1D signals with some methods (correlation, RMSE, etc) but in my opinion it limits your options in defining what it means for two waves to be similar.
EDIT
I want to clarify what phase variant/invariant means.
Here is a picture with two signals (s1 and s2), that are very similar and in phase.
I performed FFT and extracted magnitude and phase (anlge) from each Fourier coefficient. (Note, the magnitude is log compressed so the absolute values are comparable)
Magnitude for each signal (m1 and m2):
Phase for each signal (p1 and p2):
If I compute the Euclidean distance between the magnitudes and phases I get the following values:
For phase shift of 0 Magnitude distance: 6.2272 phase distance: 19.9371
I change the phase of the s2 so it will be in $\pi/2$ wrt to s1 and the values changed to:
For phase shift of 0.5 Magnitude distance: 5.9229 phase distance: 28.0612
Just for reference, here are the out of phase s1, s2, m1, m2, p1 and p2:
Signals:
Magnitude of the Fourier
Phase:
Note that the magnitude of in phase and out of phase signals are roughly unchanged. You can say that the magnitude of FFT for a signal, is insensitive to alterations of the phase (this is what the fancy "invariant" word means). The phase of FFT is another story, it is very sensitive to changes in the phase of a signal.
I post the MATLAB code so you can play with the simulations.
Hope it helps!
function FourierFeatures()
phaseShift = 0.5; % 0 == in phase, 0.5 = pi/2 out of phase
s1 = generateSignal(100,0.0);
s2 = generateSignal(100,phaseShift);
[m1 p1] = performFourier(s1);
[m2 p2] = performFourier(s2);
close all;
plot(s1,'r'); hold on; plot(s2,'b');
figure;
plot(m1,'r'); hold on; plot(m2,'b');
figure
plot(p1,'r'); hold on; plot(p2,'b');
mdist = norm(m2-m1);
pdist = norm(p2-p1);
disp(['For phase shift of ', num2str(phaseShift) , ' Magnitude distance: ',num2str(mdist), ' phase distance: ' , num2str(pdist)]);
end
%Compute FFT Magnitude and phase
function [magnitude, phase] = performFourier(signal)
f = fft(signal);
magnitude = abs(f);
magnitude = log(1+magnitude);
phase = angle(f);
end
%Generate sinusoidal signals with some noise
function signal = generateSignal(length, offset)
x = [1:length];
omega = [10, 15 ,20];
amplitude = [10 10 15];
randomAmp = [1 3 3];
y = 5 * rand(1,length);
for i = 1: size(omega,2)
y = y + amplitude(i) * sin(offset * omega(i) + x * omega(i))+ randomAmp(i) * rand(1,length);
end
signal = y;
end
