本帖最后由 highgear 于 2010-9-1 23:16 编辑
首先， 我贴出的DFT程序都是我自己写的， 而且有汇编的版本。 大家已经看到了， c 版本完全使用了16-bit 整数乘法， 32-bit 加法以及少量的移位操作， 除法（主要是 %，用于非 2的次方的 N) 可以完全避开。 可以设定在每次定时中断里采样后计算， 由于递推， 计算量很低（2个16bit乘法， 2个32bit 加法）。 唯一的问题是， 必须使用定点 scale 转换以避免浮点运算， 这不如直接使用浮点直观， 对没有处理经验的程序员可能是一个挑战。虽然演示程序为了通用起见用了 PC, 但是dft 算法程序用在 avr, 51等 8-bit mcu 是完完全全没有问题的。不要再说出单片机不能实现的胡话出来。
FFT程序我也有汇编的版本， C/C++ 版本也是采用了16-bit 整数乘法， 32-bit 加法以及少量的移位操作， 效率很高， 不过在 8-bit mcu 上可能用不上， 因为， 数据窗口点数少了， 用 dft 更好， 数据窗口点数多了， 8-bit mcu 上太慢， 不实用。 因此我就不介绍了。
最后， 我贴出我用 matlab 做的变速不变调的算法验证程序， 作为结束。
简单的讲一讲原理：
下面的程序使用了短时博立叶变换（short time fourier transform), 窗口函数为 hamming。
1） 短时博立叶变换, 这里的片断 segment = N/4, 数据被分割为 0 到 N-1, N/4 to N+N/4-1, N/2 to N+N/2-1， 依次类推。
2） 做 fft, 计算出幅度和相位。
3）计算新的幅度和相位。幅度通过插植， 相位得把 wt 计入： da(2: (1 + NX/2)) = (2*pi*segment) ./ (NX ./ (1: (NX/2)));
4）用新的幅度和相位产生新的复数， 加窗并作 ifft 生成变速后的音频数据。
SPEED = 2;
[in_rl, fs] = wavread('C:\windows\Media\Windows XP Startup.wav');
in = in_rl(:, 1)';
sizeOfData = length(in);
N = 2048;
segment = N/4;
window = hamming(N)';
X = zeros((1+ N/2), 1 + fix((sizeOfData - N)/segment));
c = 1;
for i = 0: segment: (sizeOfData - N)
fftx = fft(window .* in((i+1): (i+N)));
X(:, c) = fftx(1: (1+N/2))';
c = c + 1;
end;
[Xrows, Xcols] = size(X);
NX = 2 * (Xrows - 1);
Y = zeros(Xrows, round((Xcols - 1) / SPEED));
da = zeros(1, NX/2+1);
da(2: (1 + NX/2)) = (2*pi*segment) ./ (NX ./ (1: (NX/2)));
phase = angle(X(:, 1));
c = 1;
for i = 0: SPEED: (Xcols-2)
X1 = X(:, 1 + floor(i));
X2 = X(:, 2 + floor(i));
df = i - floor(i);
magnitude = (1-df) * abs(X1) + df * (abs(X2));
dangle = angle(X2) - angle(X1);
dangle = dangle - da' - 2 * pi * round(dangle/(2*pi));
Y(:, c) = magnitude .* exp(j * phase);
phase = phase + da' + dangle;
c = c + 1;
end
[Yrows, Ycols] = size(Y);
out = zeros(1, N + (Ycols - 1) * segment );
c = 1;
for i = 0: segment: (segment * (Ycols-1))
Yc = Y(:, c)';
Ynew = [Yc, conj(Yc((N/2): -1: 2))];
out((i+1): (i+N)) = out((i+1)i+N)) + real(ifft(Ynew)) .* window;
c = c + 1;
end;
wavplay(out, fs);
申明：以上的**和程序都是原创， 可以使用，但请不要任意转载。 | 
楼主
