二阶滤波器用C语言怎么写

2025-04-28 05:32:30
推荐回答(1个)
回答1:

这个可比你想象的复杂多了,s是个复变量,1/(s+1)极点在-1,要想用C语言写,必须理解清楚下面几个问题:
1、输入必须是个有限序列,比如(x+yi),x和y分别是两个长度为N的数组
2、要过滤的频率,必须是个整型值,或者是个整型区间
3、输出结果同样是两个长度为N的数组(p+qi)
4、整个程序需要使用最基本的复数运算,这一点C语言本身不提供,必须手工写复函数运算库
5、实现的时候具体算法还需要编,这里才是你问题的核心。

我可以送你一段FFT的程序,自己琢磨吧,和MATLAB的概念差别很大:
#include
#include
#include
#include
#include
#include
#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, 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(0.0f, 0.0f); // 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 = 0.0f, im = 0.0f; int i = 0;
for (i = 0; i < n; i++) {
rl = 0.0f; im = 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 = 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 \"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"