
1、输入必须是个有限序列,比如(x+yi),x和y分别是两个长度为N的数组
2、要过滤的频率,必须是个整型值,或者是个整型区间
3、输出结果同样是两个长度为N的数组(p+qi)
4、整个程序需要使用最基本的复数运算,这一点C语言本身不提供,必须手工写复函数运算库
5、实现的时候具体算法还需要编,这里才是你问题的核心。
我可以送你一段FFT的程序,自己琢磨吧,和MATLAB的概念差别很大:
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <windows.h>
#include "complex.h"
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, Wint n, i
long N = abs(count)long Inversing = count <0? 1: -1
for(n = 0n <N n++){ // compute from line 0 to N-1
F = CComplex(0.0f, 0.0f)// clear a line
for(i = 0i <Ni++) {
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 = 0CComplex W
for(i = 0i <N2i++) {
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 = 0i <N2i++) {
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 = 0.0f, im = 0.0fint i = 0
for (i = 0i <ni++) {
rl = 0.0fim = 0.0f
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 = 1.0f
else
factor = 1.0f / factor
//end if
for (i = 0i <ni++) {
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 \"wfft.exe\" \
(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"
float middle_filter(float middle_value [] , intcount){
float sample_value, data
int i, j
for (i=1 i for(j=count-1 j>=i,--j){
if(middle_value[j-1]=middle_value[j]{
data=middle_value[j-1]
middle_value[j-1]=middle_value[j]
middle_value[j]=data
}
}
sample_value=middle_value(count-1)/2]
return(sample_value)
}
欢迎分享,转载请注明来源:内存溢出
微信扫一扫
支付宝扫一扫
评论列表(0条)