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"
可以告诉你方法:算数平均滤波,就是求出k次采样值的总和,再除以k;中值滤波法,是把k个采样值按照从小到大排列顺序,然后找到位于最中间的那个值;最后一种差薯正不知道你们老师的防脉冲是什么意思,猜测可能是去掉k次采样值中大于虚悔或小于某个值,剩余值求平均数。手迅让别人免费给你写程序基本上不可能,这个得花时间和精力。
可以察皮的。int*** SmoothImage(int ***XImage ,int width, int height, int channel)
{
double sigma = 1.85//(n/2 -1)*0.3 +0.8 { n = 9 ,no. of elements}
double conv[3][3]
double hg = 0
//Convolution kernel
for(i=0i<3i++)
{
for(j=0j<3j++)
{
int u=i-1//含没好subtract from centre element index in 3*3 (1,1)
int v=j-1
conv[i][j] =exp ( - ((u*u) + (v*v)) / (2*sigma*sigma) )
hg += conv[i][j]
}
}
for(int i=0i<3i++)
{
for(int j=0j<3j++)
{
conv[i][j] = conv[i][j] / hg
}
}
int*** sXImage = 0
sXImage = CreateImageMatrix( sXImage , width , height , channel )// allocating a 3d array
//Assigning weights
for(int i =0i <height i++)
{
for(int j =0j <width j++)
{
for(int k =0k <谈铅 3 k++)
{
double val = 0
double valw =0
if(j-1 >0 &&i-1 >0)
{
val += conv[0][0] * XImage[i-1][j-1][k]
valw += conv[0][0]
}
if(i-1 >0 )
{
val += conv[0][1] * XImage[i-1][j][k]
valw += conv[0][1]
}
if(i-1 >0 &&j+1 <width)
{
val += conv[0][2] * XImage[i-1][j+1][k]
valw += conv[0][2]
}
if(j-1 >0 )
{
val += conv[1][0] * XImage[i][j-1][k]
valw += conv[1][0]
}
val += conv[1][1] * XImage[i][j][k]
valw += conv[1][1]
if(j+1 <width)
{
val += conv[1][2] * XImage[i][j+1][k]
valw += conv[1][2]
}
if(j-1 >0 &&i+1 <height)
{
val += conv[2][0] * XImage[i+1][j-1][k]
valw += conv[2][0]
}
if(i+1 <height)
{
val += conv[2][1] * XImage[i+1][j][k]
valw += conv[2][1]
}
if(j+1 <width &&i+1 <height)
{
val += conv[2][2] * XImage[i+1][j+1][k]
valw += conv[2][2]
}
sXImage[i][j][k] = val / valw
}
}
}
return( sXImage)
}
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)