www.pudn.com > voiceprocessingtoolbox.rar > fftViaLse01.m

% FFT via least-squares estimate 
N=8; 
if rem(N, 2)~=0 
	error('N must be even!'); 
end 
fs=1; 
time=(0:N-1)'/fs; 
x=rand(N,1)*2-1; 
 
A=ones(N,1); 
for i=1:N/2 
	A=[A, cos(2*pi*(i*fs/N)*time), sin(2*pi*(i*fs/N)*time)]; 
end 
th=A\x; 
 
plotNum=fix(N/2)+2; 
subplot(plotNum, 1, 1); 
N2=(N-1)*5+1;		% 兩點間插入四點，以便觀察波形 
timeFine=linspace(min(time), max(time), N2); 
x2=th(1)*ones(N,1); 
plot(timeFine, th(1)*ones(N2,1), '.-', time, x2, 'or');	% 畫出第一項 
ylabel('直流分量'); axis([0 N/fs -1 1]); grid on 
 
for i=1:N/2		% 畫出第二至第 1+N/2 項 
	y=th(2*i)*cos(2*pi*(i*fs/N)*time)+th(2*i+1)*sin(2*pi*(i*fs/N)*time);	% 每一項 
	x2=x2+y; 
	fprintf('i=%d, sse=%f\n', i, norm(x-x2)/sqrt(N)); 
	subplot(plotNum, 1, i+1); 
	yFine=th(2*i)*cos(2*pi*(i*fs/N)*timeFine)+th(2*i+1)*sin(2*pi*(i*fs/N)*timeFine);	% Finer verison for plotting 
	plot(timeFine, yFine, '.-', time, y, 'or'); ylabel(sprintf('諧波 %d', i)); 
	axis([0 N/fs -1 1]); grid on 
end 
 
% 畫出原來訊號 
subplot(plotNum, 1, plotNum) 
plot(time, x, 'o-'); axis([0 N/fs -1 1]); grid on 
xlabel('Time'); ylabel('原訊號'); 
 
% 將 th 轉回 fft 並比較結果 
F=fft(x); 
F2=[]; 
F2(1)=th(1)*N; 
for i=1:N/2 
	F2(i+1)=(th(2*i)-sqrt(-1)*th(2*i+1))*N/2; 
	if (i==N/2) 
		F2(i+1)=2*F2(i+1); 
	end 
end 
% fft 的對稱性 
for i=N/2+2:N 
	F2(i)=F2(N-i+2)'; 
end 
 
error1=sum(abs(F2-F.'));	% F.' 是不進行共軛轉換的轉置 
error2=sum(abs(F2-F'));		% F' 是進行共軛轉換的轉置 
fprintf('轉換成 FFT 的誤差：error1=%f, error2=%f\n', error1, error2); 
fprintf('由於 FFT 的對稱性，上述誤差應該有一項為零。\n');