DESCRIPTION OF SIGNALS
I have several sung melodies under each of nine experimental contexts ("conditions"). The melodies are all of slightly different lengths (ranging more or less from 3s to 12s), and each is encoded as a discretised melodic contour (pitch profile), at fs=1000 Hz. Those contours (rather than the original audio) represent the level of analysis of interest. See top-side of the first set of plots below.
GOALS
I'd like to characterize the average power spectrum of these melodies taken as a group (within each condition), in order to speculate about slow/fast changes in melodic contour. Some have suggested that, when applying Fourier analysis to investigate the pitch contour of a melody, the lower partials that result are responsible for the overall shape of melody, while the higher partials are responsible the more surface-level ornamentation.
Thus, it is not the frequency over time that I want (so a spectrogram per se wouldn't be the right tool), but the frequency components corresponding to each individual melodic profile, so as to later examine their average.
WHAT I TRIED
The first set of plots in the code below is a sanity check that produces the correct time-domain profiles of individual melodies (left-hand side). But their power spectrum (right-hand side) always comes up in the same shape, reminiscent of the reciprocal (1/x) function. As per the code comments, detrending the signals doesn’t help (so this is not just an unruly DC component), nor does using
FFT()
directly, instead of the Welch estimation that I chose. Is this a problem with how the spectrum was computed/displayed? Or is such a result in fact expected? Downsampling to reduce the repetition of identical values doesn't help.The second set of plots shows the average spectra shapes to ALL be (almost identical) 1/x-type shapes, across all my nine conditions. This is despite the melodic contours between and within conditions being very variable.
Although it makes sense for the averaging to be done in the frequency domain, I tried to IFFT each average back into the time domain (second set of plots, left-hand-side). The result is, once again, almost identical signals, whereas I expected to get something resembling various shapes of average melodic contours. Again, is the transform computed somehow incorrectly?
Because of the unequal melody lengths, I've had to either zero-fill up to the length of the longest melody or trim down to the length of the shortest one. None of these compromises is ideal however, as the former I believe introduces artefacts, and the latter loses me lots of samples from the longer melodies. Is there a way to keep each melody's original length?
ESSENTIAL BITS OF CODE, UNDER ONE GIVEN CONDITION
pxx_within_this_condition = [];
w_within_this_condition = [];
for i_melody = 1:N_melodies
x = retrievemel(i_melody);
x = detrend(x,0); % Detrending the input signal doesn't make a difference to the reported problem.
Nx = length(x);
window = floor(Nx * 1/2); % segment length in samples
noverlap = floor(window/2); % amount of overlap, in samples. 50% seems often used
nfft = max(256, 2^nextpow2(window));
fs = 1000;
[pxx, w] = pwelch(x, hamming(window), noverlap, nfft, 'psd'); % using just "pwelch(x)" (default), specifying spectrumtype = 'power' or 'psd', or just window instead of hamming(window) as an argument, doesnt make a difference to the reported problems
% using the plain FFT instead of the Welch estimatiom method gives artefactual-looking results also
% pxx = fft(x);
% pxx = abs(pxx/Nx);
% debug-plotting of individual melodic contours one at a time and their respective spectra
subplot(2,1,1)
plot(x)
subplot(2,1,2)
plot(pxx)
pxx_within_this_condition = [ pxx_within_this_condition; pxx' ];
w_within_this_condition = [ w_within_this_condition; w' ];
end
% average
pxx_current_condition = mean(pxx_within_this_condition, 1); % avg across the columns of the matrix (result = row vector). I understand it's ok to average only values in the magnitude/power spectrum, and ignore the phase spectrum
w_current_condition = mean( w_within_this_condition, 1);
% plot FFT as estimated by the Welch method
subplot(..)
plot(w_current_condition, 10*log10(pxx_current_condition))
% plot in the time domain the signal to which the fq-domain average above would correspond, via IFFT
subplot(..)
X = ifft(pxx_current_condition, 'symmetric'); % 'symmetric' ensures that the output is real, though not sure we really have conjugate symmetric vectors
plot(X)
PRODUCED PLOTS