一、序言
场景一：噪音信号是数据采集处理的天敌，但无时无刻它都存在，于是，信号传输时进行屏蔽防护、模数转换时给予充分的采保时间、电路实现上低通带通处理，为了减小电解电容的感抗作用有时还附加上瓷片电容滤波，有时甚至不惜损失响应时间而加大滤波系数，等等。在没有信号输入时，我们希望采集到的数据是平直的基线，但往往采集到的数据伴有严重的噪音干扰，就如下面的图形：
场景二：大部分乐器的发声频率在20HZ 到 20KHZ，但高速上行驶驾车时往往听到的不仅仅是音乐的悠扬，可能还有相当份量20KHZ以上的外部噪音。
因此需要傅氏变换获得噪音频谱，然后去掉不需要的高频噪音成份。如果对滤波后的残余噪声信号施加以相反相位的音频信号进行抵消，大致上应该是最简陋的主动降噪功效了。
二、软件实现傅氏变化
环境：Mint-Linx 21.3, freebasic 1.10, fbmath 计算库
假设信号源是200hz理想的正弦波
它受到2000hz斜纹信号干扰，相互叠加在一起，软件画出的图形是如下的样子：
对受到干扰的信号进行ADC转换并采集512个点，采集速率是2.2khz, 然后进行傅氏变换，得到的下面的频谱图形。从图上可以清楚看到2000hz处的频谱波峰，在软件中去掉1000hz以上的频谱数据，留下的就是基频信号数据了。
然后再进行反变换，将频域信号转换为时域信号：
freebasic 程序代码：
' ******************************************************************
' Fast Fourier Transform (modified from Pascal program by Don Cross)
' http://groovit.disjunkt.com/analog/time-domain/fft.html
' ******************************************************************
'
' The program generates a time signal consisting of a large 200 Hz
' sine wave added to a small 2000 Hz cosine wave, which is graphed
' on the screen (Press a key after you are done viewing each graph)
'
' Next, it performs the FFT and graphs the resulting complex
' frequency samples.
'
' Then, it filters out all frequency components above 1000 Hz in
' the transformed data.
'
' Finally, it performs the inverse transform to get a filtered
' time signal back, and graphs the result.
'
' In addition, the program compares the FFT with the direct
' computation, on a set of random data.
'
' Results are stored in the output file fftout.txt
' ******************************************************************
#INCLUDE "fbmath.bi"
#INCLUDE "fbgfx.bi"
CONST NumSamples = 512 ' Buffer size must be power of 2
CONST SamplingRate = 22050 ' Sampling rate in Hz
CONST MaxIndex = NumSamples - 1 ' Max. array index
CONST MidIndex = NumSamples \ 2
CONST DT = 1 / SamplingRate ' Time unit
CONST DF = SamplingRate / NumSamples ' Frequency unit
FUNCTION F(BYVAL T AS DOUBLE) AS DOUBLE
' Original signal
RETURN SIN(200 * TwoPi * T) + 0.2 * COS(2000 * TwoPi * T)
END FUNCTION
SUB Test_CalcFrequency(InArray() AS Complex, OutArray() AS Complex)
DIM AS Complex Z
DIM AS INTEGER I
' Fill input buffers with random data
FOR I = 0 TO MaxIndex
InArray(I).X = 10000 * RND
InArray(I).Y = 10000 * RND
NEXT I
PRINT #1, ""
PRINT #1, " *** Testing procedure CalcFrequency *** "
PRINT #1, ""
FFT InArray(), OutArray()
FOR I = 0 TO MaxIndex
Z = CalcFrequency(I, InArray())
PRINT #1, USING "##########"; I;
PRINT #1, USING "##########.######"; OutArray(I).X; Z.X; OutArray(I).Y; Z.Y
NEXT I
END SUB
SUB ListData(DataArray() AS Complex, BYREF Comment AS STRING)
DIM AS INTEGER I
PRINT #1, " *** "; Comment; " *** "
PRINT #1, ""
PRINT #1, " index real imag"
FOR I = 0 TO MaxIndex
PRINT #1, USING "##########"; I;
PRINT #1, USING "##########.######"; DataArray(I).X; DataArray(I).Y
NEXT I
PRINT #1, ""
PRINT #1, STRING(44, "-")
PRINT #1, ""
END SUB
SUB PlotData(T() AS DOUBLE, Z() AS Complex, BYREF Title AS STRING)
DIM AS DOUBLE Xmin, Xmax, Xstep ' Ox scale
DIM AS DOUBLE Ymin, Ymax, Ystep ' Oy scale
DIM AS INTEGER I ' Loop variable
DIM AS DOUBLE X(0 TO MaxIndex) ' Real part of Z
IF NOT InitScreen(12, 15, 0) THEN
PRINT "Unable to set graphic mode!"
EXIT SUB
END IF
InitGraph
FOR I = 0 TO MaxIndex : X(I) = Z(I).X : NEXT I
AutoScale T(), LinScale, Xmin, Xmax, Xstep
AutoScale X(), LinScale, Ymin, Ymax, Ystep
SetOxScale LinScale, Xmin, Xmax, Xstep
SetOyScale LinScale, Ymin, Ymax, Ystep
PlotOxAxis Ymin, 5
PlotOyAxis Xmin, 2
PlotGrid
WriteTitles Title, "Time (s)", "Amplitude"
COLOR 10 : PlotCurve T(), X(), 0, 0, -1
SLEEP
SCREEN 0
END SUB
SUB PlotFFT(Freq() AS DOUBLE, Z() AS Complex, BYREF Title AS STRING)
DIM AS DOUBLE Fq(0 TO MidIndex) ' Frequency
DIM AS DOUBLE X(0 TO MidIndex) ' Real part of FFT
DIM AS DOUBLE Y(0 TO MidIndex) ' Imag. part of FFT
' Graphic parameters for real and imaginary parts of FFT
DIM AS TCurvParam CurvParam(1 TO 2) = _
{(12, 0, 0, -1, "Real"), _
(14, 0, 0, -1, "Imag")}
DIM AS DOUBLE Xmin, Xmax, Xstep ' Ox scale
DIM AS DOUBLE Yr_min, Yr_max, Yr_step ' Oy scale (real part)
DIM AS DOUBLE Yi_min, Yi_max, Yi_step ' Oy scale (imag. part)
DIM AS DOUBLE Ymin, Ymax, Ystep ' Oy scale (global)
DIM AS INTEGER I ' Loop variable
IF NOT InitScreen(12, 15, 0) THEN
PRINT "Unable to set graphic mode!"
EXIT SUB
END IF
InitGraph
FOR I = 0 TO MidIndex
Fq(I) = Freq(I)
X(I) = Z(I).X
Y(I) = Z(I).Y
NEXT I
AutoScale Fq(), LinScale, Xmin, Xmax, Xstep
AutoScale X(), LinScale, Yr_min, Yr_max, Yr_step
AutoScale Y(), LinScale, Yi_min, Yi_max, Yi_step
#ifdef USE_PREFIX
Ymin = fbMin(Yr_min, Yi_min)
Ymax = fbMax(Yr_max, Yi_max)
Ystep = fbMin(Yr_step, Yi_step)
#else
Ymin = Min(Yr_min, Yi_min)
Ymax = Max(Yr_max, Yi_max)
Ystep = Min(Yr_step, Yi_step)
#endif
SetOxScale LinScale, Xmin, Xmax, Xstep
SetOyScale LinScale, Ymin, Ymax, Ystep
PlotOxAxis Ymin, 5
PlotOyAxis Xmin, 2
PlotGrid
WriteTitles Title, "Frequency (Hz)", "FFT"
WITH CurvParam(1)
COLOR .Col
PlotCurve Fq(), X(), .Symbol, .Size, .Connect
END WITH
WITH CurvParam(2)
COLOR .Col
PlotCurve Fq(), Y(), .Symbol, .Size, .Connect
END WITH
WriteLegend CurvParam()
SLEEP
SCREEN 0
END SUB
DIM AS DOUBLE T(0 TO MaxIndex) ' Time
DIM AS DOUBLE Freq(0 TO MaxIndex) ' Frequencies
DIM AS Complex InArray(0 TO MaxIndex) ' FFT input
DIM AS Complex OutArray(0 TO MaxIndex) ' FFT output
DIM AS INTEGER I, FreqIndex, SymIndex
FOR I = 0 TO MaxIndex
T(I) = I * DT
Freq(I) = I * DF
InArray(I).X = F(T(I))
InArray(I).Y = 0
NEXT I
OPEN "fftout.txt" FOR OUTPUT AS #1
ListData InArray(), "Time domain data before transform"
PlotData T(), InArray(), "Original signal"
FFT InArray(), OutArray()
PlotFFT Freq(), OutArray(), "Fourier Transform"
ListData OutArray(), "Frequency domain data after transform"
' Filter out everything above 1000 Hz (low-pass)
FreqIndex = 1000 / DF
SymIndex = NumSamples - FreqIndex
FOR I = 0 TO MaxIndex
IF (I > FreqIndex AND I < MidIndex) OR _
(I >= MidIndex AND I < SymIndex) THEN
OutArray(I).X = 0
OutArray(I).Y = 0
END IF
NEXT I
IFFT OutArray(), InArray()
ListData InArray(), "Time domain data after inverse transform"
PlotData T(), InArray(), "Filtered signal"
Test_CalcFrequency InArray(), OutArray()
CLOSE #1如果以上述512字节作为缓冲区，对采集的数据进行滑动(新陈代谢逐步替换)更新，可以得到连续处理后的时域信号。
注：
fbmath库是国外一个教授自己用freebasic编写的，是开源的，放在了sourceforge 上供下载。下载后用freebasic进行编译，生成libmath.a静态库，下载的包里有编写好的fbmath.bi头文件。fbmath.bi放在/usr/include/freebasic下面， fbmath.a放在/usr/lib/x86_64-linux-gnu下面(或是其它能找到的路径下面)。
fbmath包含有矩阵变换，拉氏变换，微积分方程等多种计算程序。largeint是计算超长整数的程序，比如2^65536 的计算。很好的东西，可惜平时用不上，上学时又没有。
