1.1 傅里葉變換

1.1.1 傅里葉變換概述

傅里葉變換是一種可以將滿足某個條件的函數表示成三角函數（正弦或余弦）表示的線性組合，最初，傅里葉分析是作為熱過程解析分析的工具被提出的，在電路中，這是一種分析非正弦周期信號的最便捷的工具，傅里葉是法國數學家和物理學家，在1807年的法國科學學會上發表了運用正弦曲線來描述溫度分布，但是論文中有一個具有爭議性的判定，即任何連續周期信號都可以由一組適當的正弦曲線組合而成，當時拉普拉斯與其他審查者投票通過并要發表這個論文，但是拉格朗日認為傅里葉提供的方法無法表示帶有棱角的信號，由于一些原因，這篇論文直到拉格朗日去世15年后才被發表，其實拉格朗日當時的想法是對的，正弦曲線的確無法組合一個帶有棱角的信號，但是我們可以采用正弦曲線組合成的信號來無限的逼近這個信號，直到兩種信號不存在能量差別。

通過上面的描述可以得到一個結論，傅里葉變換是一種將不規則信號通過正弦信號的線性組合表示出來的方法，但是為什么不用方波或者三角波來表示原信號呢，這是因為分解信號的方法是無窮的，但分解信號的目的是為了更加簡單地處理原來的信號。用正余弦來表示原信號會更加簡單，因為正余弦擁有原信號所不具有的性質：正弦曲線保真度。一個正弦曲線信號輸入后，輸出的仍是正弦曲線，只有幅度和相位可能發生變化，但是頻率和波的形狀仍是一樣的。且只有正弦曲線才擁有這樣的性質，正因如此我們才不用方波或三角波來表示。

我們從物理系統的特征信號角度來解釋為什么不能用其他信號的線性組合來表示原始信號，我們生活中常見的現象大多屬于 線性時不變系統 （輸入信號與輸出信號滿足線性關系，且系統參數不隨時間變換），無論你采用微分方程或者傳遞函數還是狀態空間，所以可以說正弦信號是系統的特征向量，當然指數信號也是系統的特征向量，用于標識能量的衰減與累計，系統中的衰減或擴散現象大多數是指數形式或者復指數形式，但是如果輸入的是方波或者三角波，那輸出就不一定是什么形式了，所以，除了指數信號與正弦信號外其他信號都不是線性系統的特征信號。

1.1.2 時域與頻域

從我們出生開始，我們看到的世界都是以時間貫穿，人的身高走勢，汽車的軌跡都會隨著事件發生改變，這種以時間作為參照來觀察分析動態世界的方法就屬于時域分析，頻域則是描述信號在頻率方面特性時用到的一種坐標系，頻域不是真實存在的，只是一種數學構建的領域，是一個遵循特定規則的數學范疇，正弦信號是頻域中唯一存在的波形，即正弦信號是對頻域的描述。

1.1.3 定義

也就是說我們只需要求出上面的幾個參數的值，就可以根據定義式將函數轉化為傅里葉級數的形式。與傅里葉級數相對應的就是傅里葉變換，這個定義式就簡單的多。

根據上面的兩個公式可以發現，傅里葉變換后的象函數是以頻率作為自變量的函數，而逆變換后的函數則是以時間作為自變量的函數，我們將象函數用圖像的形式表示出來就成了自變量是頻率的二維圖像，這個圖像就稱為頻譜圖。

1.2 頻域圖像

在電路分析基礎中，傅里葉經常被用于分析非正弦周期信號電路中，即講一個輸入激勵構建數學模型，采用積分的方式轉化為傅里葉級數，通過計算每個頻率狀態下的電參數最終分析電路整體的電參數，這就不僅需要最終的級數形式，還需要頻域圖像來輔助分析，在這里，我們采用MATLAB來把上面例題中的頻域圖像繪制出來，MATLAB的源代碼如下：

clc
clear
%系數計算（用于構建傅里葉級數）
syms t n;
m = 100 ;                                                    %高次諧波疊加次數
T = 2 ;                                                      %信號的周期
f = t^2 ;                                                    %函數表達式
a0 = int( f, t, -T/2, T/2 )/T ;
an = int( f*cos(n*pi*t), t, -T/2, T/2 ) ;
bn = int( f*sin(n*pi*t), t, -T/2, T/2 ) ;
%轉換定義式的形式
n = -round(T+1):1:round(T+1);
anVal = eval( an ) ;                                       %生成an系數矩陣
bnVal = eval( bn ) ;                                       %生成bn系數矩陣
An = sqrt( anVal.^2+bnVal.^2 ) ;                         %幅值開平方計算
An( round( length(n)/2 ) ) = a0 ;
phi = atan( -bnVal./anVal ) ;                             %相位計算
phi( round( length(n)/2 ) ) = 0 ;
%繪制幅度譜
figure(1);
subplot( 2, 1, 1 ) ;
stem( n, An ) ;
title( '幅度譜' ) ;
xlabel( 'n' ) ;
ylabel( 'An' ) ;
%繪制相位譜
subplot( 2, 1, 2 ) ;
plot( n, phi ) ;
title( '相位譜' ) ;
xlabel( 'n' ) ;
ylabel( '|×n' ) ;
%高次諧波擬合函數曲線
figure( 2 ) ;
t = -T: 0.01: T ;
f = a0 ;
for n=1:m
    f = f+eval(an)*cos(n*pi.*t);
    f = f+eval(bn)*sin(n*pi.*t);
end
plot( t, f ) ;
title( strcat( num2str( n ), '次諧波疊加') ) ;
xlabel( 't' ) ;
ylabel( 'f(t)' ) ;

最終仿真的結果如下圖所示。

1.3 DFT算法的MATLAB實現

傅里葉變換主要分為4類：

（1）傅里葉級數：將周期性連續函數轉換為離散頻率點

（2）連續傅里葉變換：將連續函數變換為連續頻率函數

（3）離散時間傅里葉級數DFS：將離散函數變換為連續頻率的函數

（4）離散傅里葉變換DFT：將有限長離散函數變換為離散頻率點上的函數

利用MATLAB自帶的fft函數實現DFT運算，源代碼如下。

%45點采樣
subplot( 2, 2, 1 );
N = 45 ;
n = 0:N-1 ;
t = 0.01*n ;
q = n*2*pi/N ;
x = 2*sin(4*pi*t)+5*cos(8*pi*t) ;
y = fft( x, N ) ;
plot( q, abs(y) ) ;
title('DFT N=45')
%50點采樣
subplot( 2, 2, 2 ) ;
N = 50 ;
n = 0:N-1 ;
t = 0.01*n ;
q = n*2*pi/N ;
x = 2*sin(4*pi*t)+5*cos(8*pi*t) ;
y = fft( x, N ) ;
plot( q, abs(y) ) ;
title('DFT N=50')
%55點采樣
subplot( 2, 2, 3 ) ;
N = 55 ;
n = 0:N-1 ;
t = 0.01*n ;
q = n*2*pi/N ;
x = 2*sin(4*pi*t)+5*cos(8*pi*t) ;
y = fft( x, N ) ;
plot( q, abs(y) ) ;
title('DFT N=55')
%60點采樣
subplot( 2, 2, 4 ) ;
N = 60 ;
n = 0:N-1 ;
t = 0.01*n ;
q = n*2*pi/N ;
x = 2*sin(4*pi*t)+5*cos(8*pi*t) ;
y = fft( x, N ) ;
plot( q, abs(y) ) ;
title('DFT N=60')

MATLAB運行截圖如下所示。

我們進一步增加截取長度與DFT點數，例如增加到256：

N = 256 ;
n = 0:N-1 ;
t = 0.01*n ;
q = n*2*pi/N ;
x = 2*sin(4*pi*t)+5*cos(8*pi*t) ;
y = fft( x, N ) ;
plot( q, abs(y) ) ;
title('DFT N=256')

此時信號的頻譜如下圖所示。

假設采樣頻率f=100Hz，采樣間隔T=0.01s時，取樣點N的值越大，則譜分辨率F=f/N越小，說明此時分辨率越高，加入采樣點為256時，既可以實現兩個不同頻率信號的區分。

1.4 FFT算法的MATLAB實現

FFT是離散傅立葉變換DFT的快速計算方法，適用于離散信號，并且注意變換后的點數與信號的采樣點數一致，盡管可以將信號補0，但補0不能提高頻域的分辨率。matlab中提供了函數fft做一維的FFT。

（1）時域譜和頻域譜是相互對應；時域的信號長度，決定頻域的采樣間隔，它們成導數關系；

（2）時域中信號有N點，每點間隔dt，所以時域信號長度為Ndt；那么頻譜每點的間隔就是1/(Ndt)。

clc
clear
N = 256 ;
dt = 0.02 ;
n = 0:N-1 ;
t = n*dt ;
x = sin(2*pi*t) ;
m = N ;
a = zeros( 1, m ) ;
b = zeros( 1, m ) ;
c = zeros( 1, m ) ;
for k=0:m-1
    for i=0:N-1
        a( k+1 )= a( k+1 )+2/N*x( i+1 )*cos( 2*pi*k*i/N ) ;
        b( k+1 )= b( k+1 )+2/N*x( i+1 )*sin( 2*pi*k*i/N ) ;
    end
c( k+1 ) = sqrt( a( k+1 )^2+b( k+1 )^2 ) ;
end
%時域圖
subplot( 2, 1, 1 ) ;
plot( t, x ) ;
title( '原始信號' ) ;
xlabel( 't' ) ;
ylabel( 'f(t)' ) ;
%頻域圖像
f = ( 0:m-1 )/( N*dt ) ;
subplot( 2, 1, 2 ) ;
plot( f, c ) ;hold on
title( 'Fourier' ) ;
xlabel( '頻率/HZ' ) ;
ylabel( '幅值' ) ;
%計算頻率
ind = find( c==max(c), 1, 'first' ) ;                                  %查找頻率第一個最大值
x0 = f( ind );                                                              %得到橫坐標
y0 = c( ind );                                                              %得到縱坐標
plot( x0, y0, 'ro' );hold off
text( x0+1, y0-0.1, num2str( x0, '頻率=%f' ) ) ;

MATLAB運行結果如下圖所示。

1.3基于51單片機的頻譜設計

1.3.1 功能

我們現在使用51單片機來實現頻譜的顯示，通過輸入音頻信號來輸出對應的頻段的幅值，這次選用的是STC12C5A60S2，運行時鐘33MHz，內置ADC采樣，利用LED組成點陣顯示屏幕，通過ADC不停的采樣音頻數據，將數據轉換為頻域幅值顯示在LED上面。

1.3.2 原理圖分析

（1）最小系統電路

單片機系統采用國產的51單片機，在之前的51單片機中我們已經具體了解過51單片機的應用，這里不再詳細贅述，我們這次采用的是33MHz的晶振，STC12C系列相對于傳統的51單片機來說，運行速度快了接近12倍，可以運行256位的FFT算法，程序采用串口的方式燒錄進去，其中P1口可以設置為ADC輸入口，內置8通道10位ADC轉換器，采用直流5V供電。

（2）LED驅動電路

由于我們使用的單片機只有44個引腳，實際能用的IO端口僅有35個，而我們是用的是17×16的矩陣，需要33個端口，去掉1組串口，1個ADC，剩下的端口速度各不相同，已經無法滿足驅動的需要，所以采用74HC595芯片作為驅動電路，74HC595是一個可以將串行數據轉換為并行數據的芯片，并且可以通過級聯的方式輸出數據，我們采用2片595芯片的級聯來實現串行數據轉16位并行的的效果，其余的1位數據通過IO口來達到控制的目的。74HC595的時序圖如下圖所示。

（3）點陣顯示電路

這里采用17×16的點陣作為顯示，這種點陣連接的方式可以通過控制17+16=33個端口來實現17*16=272個LED亮滅的控制。