freq.c

MIDI解码程序(用VC编写)

    }

    /* go out 2 zero crossings in both directions, starting at maxpos */
    /* find the peaks after the 2nd crossing */
    minoffset1 = 0;
    for (n = 0, oldamp = origdata[maxpos], i = maxpos - 1; i >= 0 && n < 2; i--)
    {
	amp = origdata[i];
	if ((oldamp && amp == 0) || (oldamp > 0 && amp < 0) ||
	    (oldamp < 0 && amp > 0))
	    n++;
	oldamp = amp;
    }
    minoffset1 = i;
    maxamp = labs(origdata[i]);
    while (i >= 0)
    {
	amp = origdata[i];
	if ((oldamp && amp == 0) || (oldamp > 0 && amp < 0) ||
	    (oldamp < 0 && amp > 0))
	{
	    break;
	}
	oldamp = amp;

	amp = labs(amp);
	if (amp > maxamp)
	{
	    maxamp = amp;
	    minoffset1 = i;
	}
	i--;
    }

    minoffset2 = 0;
    for (n = 0, oldamp = origdata[maxpos], i = maxpos + 1; i < length0 && n < 2; i++)
    {
	amp = origdata[i];
	if ((oldamp && amp == 0) || (oldamp > 0 && amp < 0) ||
	    (oldamp < 0 && amp > 0))
	    n++;
	oldamp = amp;
    }
    minoffset2 = i;
    maxamp = labs(origdata[i]);
    while (i < length0)
    {
	amp = origdata[i];
	if ((oldamp && amp == 0) || (oldamp > 0 && amp < 0) ||
	    (oldamp < 0 && amp > 0))
	{
	    break;
	}
	oldamp = amp;

	amp = labs(amp);
	if (amp > maxamp)
	{
	    maxamp = amp;
	    minoffset2 = i;
	}
	i++;
    }

    /* upper bound on the detected frequency */
    /* distance between the two peaks is at most 2 periods */
    minoffset = (minoffset2 - minoffset1);
    if (minoffset < 4)
	minoffset = 4;
    max_guessfreq = (float) rate / (minoffset * 0.5);
    if (max_guessfreq >= (rate >> 1)) max_guessfreq = (rate >> 1) - 1;

    /* lower bound on the detected frequency */
    maxoffset = rate / pitch_freq_lb_table[LOWEST_PITCH] + 0.5;
    if (maxoffset > (length >> 1))
	maxoffset = (length >> 1);
    min_guessfreq = (float) rate / maxoffset;

    /* perform the in place FFT */
    rdft(length, 1, floatdata, ip, w);

    /* calc the magnitudes */
    for (i = 2; i < length; i++)
    {
	mag = floatdata[i++];
	mag *= mag;
	mag += floatdata[i] * floatdata[i];
	magdata[i >> 1] = sqrt(mag);
    }

    /* find max mag */
    maxmag = 0;
    for (i = 1; i < (length >> 1); i++)
    {
	mag = magdata[i];

	pitch = fft1_bin_to_pitch[i];
	if (pitch && mag > maxmag)
	    maxmag = mag;
    }
    
    /* Apply non-linear scaling to the magnitudes
     * I don't know why this improves the pitch detection, but it does
     * The best choice of power seems to be between 1.64 - 1.68
     */
    for (i = 1; i < (length >> 1); i++)
	magdata[i] = maxmag * pow(magdata[i] / maxmag, 1.66);

    /* bin the pitches */
    for (i = 1; i < (length >> 1); i++)
    {
	mag = magdata[i];

	pitch = fft1_bin_to_pitch[i];
	pitchbins[pitch] += mag;

	if (mag > pitchmags[pitch])
	    pitchmags[pitch] = mag;
    }

    /* zero out lowest pitch, since it contains all lower frequencies too */
    pitchbins[LOWEST_PITCH] = 0;

    /* find the largest peak */
    for (i = LOWEST_PITCH + 1, maxsum = -42; i <= HIGHEST_PITCH; i++)
    {
	sum = pitchbins[i];
	if (sum > maxsum)
	{
	    maxsum = sum;
	    largest_peak = i;
	}
    }

    minpitch = assign_pitch_to_freq(min_guessfreq);
    if (minpitch > HIGHEST_PITCH) minpitch = HIGHEST_PITCH;

    /* zero out any peak below minpitch */
    for (i = LOWEST_PITCH + 1; i < minpitch; i++)
	pitchbins[i] = 0;

    /* remove all pitches below threshold */
    for (i = minpitch; i <= HIGHEST_PITCH; i++)
    {
	if (pitchbins[i] / maxsum < 0.01 && pitchmags[i] / maxmag < 0.01)
	    pitchbins[i] = 0;
    }

    /* keep local maxima */
    for (i = LOWEST_PITCH + 1; i < HIGHEST_PITCH; i++)
    {
    	double temp;

	temp = pitchbins[i];
	    
    	/* also keep significant bands to either side */
    	if (temp && pitchbins[i-1] < temp && pitchbins[i+1] < temp)
    	{
    	    new_pitchbins[i] = temp;

    	    temp *= 0.5;
    	    if (pitchbins[i-1] >= temp)
    	    	new_pitchbins[i-1] = pitchbins[i-1];
    	    if (pitchbins[i+1] >= temp)
    	    	new_pitchbins[i+1] = pitchbins[i-1];
    	}
    }
    memcpy(pitchbins, new_pitchbins, 129 * sizeof(double));

    /* find lowest and highest pitches */
    minpitch = LOWEST_PITCH;
    while (minpitch < HIGHEST_PITCH && !pitchbins[minpitch])
    	minpitch++;
    maxpitch = HIGHEST_PITCH;
    while (maxpitch > LOWEST_PITCH && !pitchbins[maxpitch])
    	maxpitch--;

    /* uh oh, no pitches left...
     * best guess is middle C
     * return 260 Hz, since exactly 260 Hz is never returned except on error
     * this should only occur on blank/silent samples
     */
    if (maxpitch < minpitch)
    {
    	free(floatdata);
    	return 260.0;
    }

    /* pitch assignment bounds based on zero crossings and pitches kept */
    if (pitch_freq_lb_table[minpitch] > min_guessfreq)
    	min_guessfreq = pitch_freq_lb_table[minpitch];
    if (pitch_freq_ub_table[maxpitch] < max_guessfreq)
    	max_guessfreq = pitch_freq_ub_table[maxpitch];

    minfreq = pitch_freq_lb_table[minpitch];
    if (minfreq >= (rate >> 1)) minfreq = (rate >> 1) - 1;
    
    maxfreq = pitch_freq_ub_table[maxpitch];
    if (maxfreq >= (rate >> 1)) maxfreq = (rate >> 1) - 1;

    minbin = minfreq / f0;
    if (!minbin)
	minbin = 1;
    maxbin = ceil(maxfreq / f0);
    if (maxbin >= (length >> 1))
    	maxbin = (length >> 1) - 1;

    /* filter out all "noise" from magnitude array */
    for (i = minbin, n = 0; i <= maxbin; i++)
    {
	pitch = fft1_bin_to_pitch[i];
	if (pitchbins[pitch])
	{
	    prunemagdata[i] = magdata[i];
	    n++;
	}
    }

    /* whoa!, there aren't any strong peaks at all !!! bomb early
     * best guess is middle C
     * return 260 Hz, since exactly 260 Hz is never returned except on error
     * this should only occur on blank/silent samples
     */
    if (!n)
    {
    	free(floatdata);
	return 260.0;
    }

    memset(new_pitchbins, 0, 129 * sizeof(double));

    maxsum = -1;
    minpitch = assign_pitch_to_freq(min_guessfreq);
    maxpitch = assign_pitch_to_freq(max_guessfreq);
    maxpitch2 = assign_pitch_to_freq(max_guessfreq) + 9;
    if (maxpitch2 > HIGHEST_PITCH) maxpitch2 = HIGHEST_PITCH;

    /* initial guess is first local maximum */
    bestfreq = pitch_freq_table[minpitch];
    if (minpitch < HIGHEST_PITCH &&
    	pitchbins[minpitch+1] > pitchbins[minpitch])
    	    bestfreq = pitch_freq_table[minpitch+1];

    /* find best fundamental */
    for (i = minpitch; i <= maxpitch2; i++)
    {
	if (!pitchbins[i])
	    continue;

    	minfreq2 = pitch_freq_lb_table[i];
    	maxfreq2 = pitch_freq_ub_table[i];
    	freq_inc = (maxfreq2 - minfreq2) * 0.1;
    	if (minfreq2 >= (rate >> 1)) minfreq2 = (rate >> 1) - 1;
    	if (maxfreq2 >= (rate >> 1)) maxfreq2 = (rate >> 1) - 1;

	/* look for harmonics */
	for (freq = minfreq2; freq <= maxfreq2; freq += freq_inc)
	{
	    double ratio;
	
    	    n = total = 0;
    	    sum = weightsum = 0;

	    for (j = 1; j <= 32 && (newfreq = j*freq) <= maxfreq; j++)
    	    {
    	    	pitch = assign_pitch_to_freq(newfreq);

    	    	if (pitchbins[pitch])
    	    	{
    	    	    sum += pitchbins[pitch];
    	    	    n++;
    	    	    total = j;
    	  	}
    	    }

	    /* only pitches with good harmonics are assignment candidates */
	    if (n > 1)
	    {
	    	double ratio;
	    	
	    	ratio = (double) n / total;
	    	if (ratio >= 0.333333)
	    	{
	    	    weightsum = ratio * sum;
	    	    pitch = assign_pitch_to_freq(freq);

		    /* use only these pitches for chord detection */
	    	    if (pitch <= HIGHEST_PITCH && pitchbins[pitch])
	    	    	new_pitchbins[pitch] = weightsum;

		    if (pitch > maxpitch)
		    	continue;

    	    	    if (n < 2 || weightsum > maxsum)
    	    	    {
    	    	    	maxsum = weightsum;
    	    	    	bestfreq = freq;
    	    	    }
    	    	}
    	    }
    	}
    }

    bestpitch = assign_pitch_to_freq(bestfreq);

    /* assign chords */
    if ((pitch = assign_chord(new_pitchbins, chord,
    	bestpitch - 9, maxpitch2, bestpitch)) >= 0)
    	    bestpitch = pitch;

    bestfreq = pitch_freq_table[bestpitch];

    /* tune based on the fundamental and harmonics up to +5 octaves */
    sum = weightsum = 0;
    for (i = 1; i <= 32 && (freq = i*bestfreq) <= maxfreq; i++)
    {
	double tune;

	minfreq2 = pitch_freq_lb_table[bestpitch];
	maxfreq2 = pitch_freq_ub_table[bestpitch];

	minbin = minfreq2 * f0_inv;
	if (!minbin) minbin = 1;
	maxbin = ceil(maxfreq2 * f0_inv);
	if (maxbin >= (length>>1))
	    maxbin = (length>>1) - 1;

	for (bin = minbin; bin <= maxbin; bin++)
	{
	    tune = -36.37631656 + 17.31234049 * log(bin*f0) - bestpitch;
	    sum += magdata[bin];
	    weightsum += magdata[bin] * tune;
	}
    }

    bestfreq = 13.75 * exp(((bestpitch + weightsum / sum) - 9) /
    	       12 * log(2));

    /* Since we are using exactly 260 Hz as an error code, fudge the freq
     * on the extremely unlikely chance that the detected pitch is exactly
     * 260 Hz.
     */
    if (bestfreq == 260.0)
    	bestfreq += 1E-5;

    free(floatdata);

    return bestfreq;
}



int assign_pitch_to_freq(float freq)
{

    int pitch;

    /* round to nearest integer using: ceil(fraction - 0.5) */
    /* -0.5 is already added into the first constant below */
    pitch = ceil(-36.87631656f + 17.31234049f * log(freq));

    /* min and max pitches */
    if (pitch < LOWEST_PITCH) pitch = LOWEST_PITCH;
    else if (pitch > HIGHEST_PITCH) pitch = HIGHEST_PITCH;

    return pitch;
}