Artefactual results in the Welch PSD (FFT) of melodic contours (pitch profiles)

Question

DESCRIPTION OF SIGNALS

I have several sung melodies under each of nine 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 one sample per ms (fs=1000 Hz). It is those, rather than the original audio, that represent the level of analysis of interest. See top-side of the first set of plots below.

GOALS

First, 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, such as Marc Evans - Spectral Melody Decomposition 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.

Second, I'd like to compare these characterizations across conditions.

PROBLEMS ENCOUNTERED & FIXES SOUGHT

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? I tried downsampling to reduce the repetition of identical values, but it doesn't make a difference. I also thought that the signal consisting of only a few (very slow) step changes, rather than a proper oscillation of suitable length, might be a problem; but computing the spectrum for a large contour (obtained by merging together all contours within one condition) produces a similar result.

2. Worse still, 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 (COMMENTED, MATLAB) 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