matlab中fft的使用问题

貌似这是个书呆子编滴程序，以为乘法有交换律就随意改变位置，完全不管物理意义好不好理解，其实不就是采样定理要除2嘛！fft后最高频率是采样频率滴一半，计算频率数组时是以此原理再按几分之几均分其它部分，其实最易理解滴是(0:length(Yc)-1) /length(Yc)(Fs/2)这种表达，只不过原式与这式数学上计算结果相同。还有乘除法不用点运算，小心编译出错，真是相当糟糕的代码！

这个可比你想象的复杂多了，s是个复变量，1/(s+1)极点在-1，要想用C语言写，必须理解清楚下面几个问题：
1、输入必须是个有限序列，比如(x+yi)，x和y分别是两个长度为N的数组
2、要过滤的频率，必须是个整型值，或者是个整型区间
3、输出结果同样是两个长度为N的数组(p+qi)
4、整个程序需要使用最基本的复数运算，这一点C语言本身不提供，必须手工写复函数运算库
5、实现的时候具体算法还需要编，这里才是你问题的核心。
我可以送你一段FFT的程序，自己琢磨吧，和MATLAB的概念差别很大：
#include <asserth>
#include <mathh>
#include <stdioh>
#include <stdlibh>
#include <stringh>
#include <windowsh>
#include "complexh"
extern "C" {
// Discrete Fourier Transform (Basic Version, Without Any Enhancement)
// return - Without Special Meaning, constantly, zero
int DFT (long count, CComplex input, CComplex output)
{
assert(count);
assert(input);
assert(output);
CComplex F, X, T, W; int n, i;
long N = abs(count); long Inversing = count < 0 1: -1;
for(n = 0; n < N ; n++){ // compute from line 0 to N-1
F = CComplex(00f, 00f); // clear a line
for(i = 0; i < N; i++) {
T = input[i];

W = HarmonicPI2(Inversing n i, N);
X = T W;
F += X; // fininshing a line
}//next i
// save data to outpus
memcpy(output + n, &F, sizeof(F));
}//next n
return 0;
}//end DFT
int fft (long count, CComplex input, CComplex output)
{
assert(count);
assert(input);
assert(output);
int N = abs(count); long Inversing = count < 0 -1: 1;
if (N % 2 || N < 5) return DFT(count, input, output);
long N2 = N / 2;
CComplex iEven = new CComplex[N2]; memset(iEven, 0, sizeof(CComplex) N2);
CComplex oEven = new CComplex[N2]; memset(oEven, 0, sizeof(CComplex) N2);
CComplex iOdd = new CComplex[N2]; memset(iOdd , 0, sizeof(CComplex) N2);
CComplex oOdd = new CComplex[N2]; memset(oOdd , 0, sizeof(CComplex) N2);
int i = 0; CComplex W;
for(i = 0; i < N2; i++) {
iEven[i] = input[i 2];
iOdd [i] = input[i 2 + 1];
}//next i
fft(N2 Inversing, iEven, oEven);
fft(N2 Inversing, iOdd, oOdd );
for(i = 0; i < N2; i++) {
W = HarmonicPI2(Inversing (- i), N);
output[i] = oEven[i] + W oOdd[i];
output[i + N2] = oEven[i] - W oOdd[i];
}//next i
return 0;
}//end FFT
void __stdcall FFT(
long N, // Serial Length, N > 0 for DFT, N < 0 for iDFT - inversed Discrete Fourier Transform
double inputReal, double inputImaginary, // inputs
double AmplitudeFrequences, double PhaseFrequences) // outputs
{
if (N == 0) return;
if (!inputReal && !inputImaginary) return;
short n = abs(N);
CComplex input = new CComplex[n]; memset(input, 0, sizeof(CComplex) n);
CComplex output= new CComplex[n]; memset(output,0, sizeof(CComplex) n);
double rl = 00f, im = 00f; int i = 0;
for (i = 0; i < n; i++) {
rl = 00f; im = 00f;
if (inputReal) rl = inputReal[i];
if (inputImaginary) im = inputImaginary[i];
input[i] = CComplex(rl, im);
}//next i
int f = fft(N, input, output);
double factor = n;
//factor = sqrt(factor);
if (N > 0)
factor = 10f;
else
factor = 10f / factor;
//end if
for (i = 0; i < n; i++) {
if (AmplitudeFrequences) AmplitudeFrequences[i] = output[i]getReal() factor;
if (PhaseFrequences) PhaseFrequences[i] = output[i]getImaginary() factor;
}//next i
delete [] output;
delete [] input;
return ;
}//end FFT
int __cdecl main(int argc, char argv[])
{
fprintf(stderr, "%s usage:\n", argv[0]);
fprintf(stderr, "Public Declare Sub FFT Lib \"wfftexe\" \
(ByVal N As Long, ByRef inputReal As Double, ByRef inputImaginary As Double, \
ByRef freqAmplitude As Double, ByRef freqPhase As Double)");
return 0;
}//end main
};//end extern "C"

t=0:1/256:1;%采样步长
y= 2+3cos(2pi50t-pi30/180)+15cos(2pi75t+pi90/180);
N=length(t); %样点个数
plot(t,y);
fs=256;%采样频率
df=fs/(N-1);%分辨率
f=(0:N-1)df;%其中每点的频率
Y=fft(y(1:N))/N2;%真实的幅值
%Y=fftshift(Y);
figure(2)
plot(f(1:N/2),abs(Y(1:N/2)));

IDFT就是Inverse Discrete Fourier Transform 离散傅里叶逆变换。
FFT就是Fast Fourier Transform 快速傅里叶变换。
两者的应用都是将时域中难以处理的信号转换成易于复处理的频域信号，分析完成后进行傅里叶反变换即得到原始的时域信号。
两者的异同是：我们知道在数学上用级数来无限逼进某个函数，以便简化计算过程而又不致使误差过大，这样工程上才能应用，否则一些数学模型是无法实现快速求解的。
IDFT：对于有限长的序列我们可以使用离散傅立叶变换，IDFT是对制序列傅立叶变换的等距采样。zhidao
FFT：并不是与IDFT不相同的另一种变换（即原理是一样的），而是为了减少IDFT运算次数的一种快速算法。它是对IDFT变换式进行一次次的分解，使其成为若干小点数IDFT的组合，从而减小运算量。常用的FFT是以2为基数，它的运算效率高，程序比较简单，使用也十分地方便。

姓名：王柯祎

学号：20021110373T

转自 ：>