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.


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.


  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?


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
    pxx_within_this_condition = [   pxx_within_this_condition;      pxx'    ];
    w_within_this_condition = [ w_within_this_condition;            w'  ];  

% 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
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
X = ifft(pxx_current_condition, 'symmetric'); % 'symmetric' ensures that the output is real, though not sure we really have conjugate symmetric vectors


enter image description here

enter image description here

enter image description here

  • $\begingroup$ Sorry, this is very difficult to understand. Why would you take the spectrum of the pitch profile? The pitch profile consists mainly of flat lines, so low frequencies will dominate by a wide margin. What did you expect to see? What exactly is the metric you are looking for. What units does it have and what should it be for the example that you have plotted. It could be "number of pitch changes per second", "weighted average of pitch changes per second". etc. Can you clarify? $\endgroup$
    – Hilmar
    Apr 21, 2023 at 19:08
  • $\begingroup$ Apologies, I realise the question is rather complex and unclear. Why take the spectrum of the pitch profile and what do I expect to see: the article I linked under the Goals section also deals with 'plateaus' (melodic contour made up of discrete pitches), and found differences among lower and higher partials that result – so it's those kind of differences I'm hoping to find. "Rate (or: patterns) of pitch change in a melody" (largely construed) could be the metric, but I wanted this to be a rather botom-up exploration. $\endgroup$
    – z8080
    Apr 22, 2023 at 12:02
  • $\begingroup$ Also, the plot example I gave of one individual melodic profile wasn't the best as note durations were - in a special case - all the same, so I've updated this with two more representative profiles. I expected their PSD to differ somehow, but PSDs always looks like noisy 1/x functions. $\endgroup$
    – z8080
    Apr 22, 2023 at 12:02

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.