Power spectral density of melodic contours

Question

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

1. 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.

2. 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.

3. 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?

4. 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

Answer 1

The $1/f$ behavior is fairly easy to explain. Let's consider a single pitch profile $x[n]$ and it's first derivative $x'[n]$.

Since $x[n]$ is a step function that's mostly flat, $x'[n]$ is zero except at the step locations so we can express it as

$$x'[n] = \sum_{k=0}^{K-1} s_k\cdot\delta[n-n_k]$$

where K is the number of steps, $s_k$ the height of step $k$ and $n_k$ the location in time of step $k$. $\delta[n]$ is the unit impulse response, the PSD of which is white, i.e. $DFT\{\delta[n]\} = 1$.

Unless the individual steps are highly correlated both in time and in amplitude, the sum of the unit impulses will also have a white spectrum (with a few noise-like wiggles on top).

To get back from the derivative to the original signal we need to integrate, which corresponds to multiplication with $1/j\omega$ in the frequency domain. Hence the spectrum of these step function will always have a $1/f$ shape with a few small wiggles. The only difference that the exact profile makes will be the fine structure of the wiggles, but these are small and will be hard to interpret.

I don't fully understand what you are trying to measure, but I don't the think PSD of the pitch profile is the way to go here. You can try looking at the PSD of the first derivative of the pitch profile: this way you just see the wiggle without the $1/f$ on top. I have no idea whether this is useful or not.

In order to determine the best algorithm you need to first clearly define the metric. For example if you want the amount of "pitch variability" you could simply take the standard deviation of the pitch profile. You can also calculate the step frequency (how many pitch changes per second) the average pitch step size, etc. But without knowing exactly what you want, its hard to tell which one is the best choice.