DSP 28335 傅里叶变换和滤波程序实现
1. 系统配置和头文件
#include "DSP28x_Project.h" // DSP2833x头文件
#include "math.h"
#include "fpu.h"
// 定义FFT相关参数
#define FFT_SIZE 1024 // FFT点数
#define LOG2_FFT_SIZE 10 // log2(FFT_SIZE)
#define PI 3.14159265358979
// 全局变量声明
#pragma DATA_SECTION(fft_input, "FFTBuffer");
#pragma DATA_SECTION(fft_output, "FFTBuffer");
#pragma DATA_SECTION(fir_coeffs, "FIRCoef");
#pragma DATA_SECTION(iir_coeffs, "IIRCoef");
float fft_input[FFT_SIZE]; // FFT输入缓冲区
float fft_output[FFT_SIZE]; // FFT输出缓冲区
float fft_mag[FFT_SIZE/2]; // 幅值谱
float fft_phase[FFT_SIZE/2]; // 相位谱
// 窗函数
float window[FFT_SIZE];
2. FFT实现（基2时间抽取）
// 复数结构体
typedef struct {
float real;
float imag;
} Complex;
// 位反转函数
void bit_reverse(Complex* data, int n)
{
int i, j, k;
Complex temp;
j = 0;
for(i = 0; i < n-1; i++)
{
if(i < j)
{
temp = data[i];
data[i] = data[j];
data[j] = temp;
}
k = n >> 1;
while(k <= j)
{
j -= k;
k >>= 1;
}
j += k;
}
}
// 复数乘法
Complex complex_multiply(Complex a, Complex b)
{
Complex result;
result.real = a.real * b.real - a.imag * b.imag;
result.imag = a.real * b.imag + a.imag * b.real;
return result;
}
// 蝶形运算
void butterfly(Complex* data, int n, int m)
{
int i, j, k, le, le2;
Complex u, w, t;
le = 1 << m;
le2 = le >> 1;
// 初始化旋转因子
u.real = 1.0;
u.imag = 0.0;
w.real = cos(PI / le2);
w.imag = -sin(PI / le2);
for(j = 0; j < le2; j++)
{
for(i = j; i < n; i += le)
{
int ip = i + le2;
t = complex_multiply(data[ip], u);
data[ip].real = data[i].real - t.real;
data[ip].imag = data[i].imag - t.imag;
data[i].real = data[i].real + t.real;
data[i].imag = data[i].imag + t.imag;
}
u = complex_multiply(u, w);
}
}
// 主FFT函数
void fft(float* input, float* output_real, float* output_imag, int n, int log2n)
{
int i;
Complex* data = (Complex*)output_real; // 复用输出缓冲区
// 准备复数数据
for(i = 0; i < n; i++)
{
data[i].real = input[i];
data[i].imag = 0.0;
}
// 位反转
bit_reverse(data, n);
// 蝶形运算
for(i = 1; i <= log2n; i++)
{
butterfly(data, n, i);
}
// 分离实部和虚部
for(i = 0; i < n; i++)
{
output_real[i] = data[i].real;
output_imag[i] = data[i].imag;
}
}
3. 窗函数生成
// 汉宁窗
void hanning_window(float* window, int n)
{
int i;
for(i = 0; i < n; i++)
{
window[i] = 0.5 * (1.0 - cos(2.0 * PI * i / (n - 1)));
}
}
// 海明窗
void hamming_window(float* window, int n)
{
int i;
for(i = 0; i < n; i++)
{
window[i] = 0.54 - 0.46 * cos(2.0 * PI * i / (n - 1));
}
}
// 应用窗函数
void apply_window(float* data, float* window, int n)
{
int i;
for(i = 0; i < n; i++)
{
data[i] *= window[i];
}
}
4. 幅值和相位计算
// 计算幅值谱
void compute_magnitude(float* real, float* imag, float* mag, int n)
{
int i;
for(i = 0; i < n/2; i++) // 只计算前一半（对称）
{
mag[i] = sqrt(real[i] * real[i] + imag[i] * imag[i]);
}
}
// 计算相位谱
void compute_phase(float* real, float* imag, float* phase, int n)
{
int i;
for(i = 0; i < n/2; i++)
{
if(real[i] != 0.0)
{
phase[i] = atan2(imag[i], real[i]) * 180.0 / PI;
}
else
{
phase[i] = 0.0;
}
}
}
5. FIR滤波器实现
// FIR滤波器结构体
typedef struct {
float* coeffs; // 滤波器系数
float* buffer; // 延迟线缓冲区
int length; // 滤波器长度
int index; // 当前索引
} FIRFilter;
// 初始化FIR滤波器
void fir_init(FIRFilter* fir, float* coeffs, float* buffer, int length)
{
fir->coeffs = coeffs;
fir->buffer = buffer;
fir->length = length;
fir->index = 0;
// 清空缓冲区
int i;
for(i = 0; i < length; i++)
{
fir->buffer[i] = 0.0;
}
}
// FIR滤波器处理
float fir_filter(FIRFilter* fir, float input)
{
int i;
float output = 0.0;
// 更新缓冲区
fir->buffer[fir->index] = input;
// 计算卷积
for(i = 0; i < fir->length; i++)
{
int idx = (fir->index - i + fir->length) % fir->length;
output += fir->coeffs[i] * fir->buffer[idx];
}
// 更新索引
fir->index = (fir->index + 1) % fir->length;
return output;
}
// 批量FIR滤波
void fir_filter_batch(FIRFilter* fir, float* input, float* output, int n)
{
int i;
for(i = 0; i < n; i++)
{
output[i] = fir_filter(fir, input[i]);
}
}
6. IIR滤波器实现
// 二阶IIR滤波器结构体
typedef struct {
float a1, a2; // 分母系数
float b0, b1, b2; // 分子系数
float x1, x2; // 输入延迟
float y1, y2; // 输出延迟
} IIRBiquad;
// 初始化二阶IIR滤波器
void iir_biquad_init(IIRBiquad* iir, float b0, float b1, float b2, float a1, float a2)
{
iir->b0 = b0;
iir->b1 = b1;
iir->b2 = b2;
iir->a1 = a1;
iir->a2 = a2;
iir->x1 = 0.0;
iir->x2 = 0.0;
iir->y1 = 0.0;
iir->y2 = 0.0;
}
// 二阶IIR滤波器处理
float iir_biquad_filter(IIRBiquad* iir, float input)
{
float output = iir->b0 * input + iir->b1 * iir->x1 + iir->b2 * iir->x2
- iir->a1 * iir->y1 - iir->a2 * iir->y2;
// 更新延迟
iir->x2 = iir->x1;
iir->x1 = input;
iir->y2 = iir->y1;
iir->y1 = output;
return output;
}
7. 滤波器设计函数
// 生成低通FIR滤波器系数（窗函数法）
void design_fir_lpf(float* coeffs, int length, float fc, float fs)
{
int i, center;
float omega_c = 2.0 * PI * fc / fs;
center = length / 2;
for(i = 0; i < length; i++)
{
if(i == center)
{
coeffs[i] = omega_c / PI;
}
else
{
float n = i - center;
coeffs[i] = sin(omega_c * n) / (PI * n);
}
// 应用汉宁窗
coeffs[i] *= 0.5 * (1.0 - cos(2.0 * PI * i / (length - 1)));
}
}
// 生成带通FIR滤波器系数
void design_fir_bpf(float* coeffs, int length, float f1, float f2, float fs)
{
int i, center;
float omega1 = 2.0 * PI * f1 / fs;
float omega2 = 2.0 * PI * f2 / fs;
center = length / 2;
for(i = 0; i < length; i++)
{
if(i == center)
{
coeffs[i] = (omega2 - omega1) / PI;
}
else
{
float n = i - center;
coeffs[i] = (sin(omega2 * n) - sin(omega1 * n)) / (PI * n);
}
// 应用海明窗
coeffs[i] *= 0.54 - 0.46 * cos(2.0 * PI * i / (length - 1));
}
}
8. 实时处理模块
// 实时FFT处理结构体
typedef struct {
float input_buffer[FFT_SIZE];
float fft_real[FFT_SIZE];
float fft_imag[FFT_SIZE];
float mag[FFT_SIZE/2];
float phase[FFT_SIZE/2];
int buffer_index;
int ready;
} RealTimeFFT;
// 初始化实时FFT
void realtime_fft_init(RealTimeFFT* rt_fft)
{
int i;
for(i = 0; i < FFT_SIZE; i++)
{
rt_fft->input_buffer[i] = 0.0;
rt_fft->fft_real[i] = 0.0;
rt_fft->fft_imag[i] = 0.0;
}
for(i = 0; i < FFT_SIZE/2; i++)
{
rt_fft->mag[i] = 0.0;
rt_fft->phase[i] = 0.0;
}
rt_fft->buffer_index = 0;
rt_fft->ready = 0;
}
// 添加新样本
void realtime_fft_add_sample(RealTimeFFT* rt_fft, float sample)
{
rt_fft->input_buffer[rt_fft->buffer_index] = sample;
rt_fft->buffer_index++;
if(rt_fft->buffer_index >= FFT_SIZE)
{
rt_fft->buffer_index = 0;
rt_fft->ready = 1; // 缓冲区已满，可以执行FFT
}
}
// 执行实时FFT
void realtime_fft_process(RealTimeFFT* rt_fft)
{
if(rt_fft->ready)
{
// 应用窗函数
apply_window(rt_fft->input_buffer, window, FFT_SIZE);
// 执行FFT
fft(rt_fft->input_buffer, rt_fft->fft_real, rt_fft->fft_imag, FFT_SIZE, LOG2_FFT_SIZE);
// 计算幅值和相位
compute_magnitude(rt_fft->fft_real, rt_fft->fft_imag, rt_fft->mag, FFT_SIZE);
compute_phase(rt_fft->fft_real, rt_fft->fft_imag, rt_fft->phase, FFT_SIZE);
rt_fft->ready = 0;
}
}
9. 主程序和ADC配置
// 全局变量
FIRFilter fir_lpf;
IIRBiquad iir_biquad;
RealTimeFFT rt_fft;
float fir_coeffs[65]; // 65阶FIR滤波器系数
float fir_buffer[65]; // FIR滤波器延迟线
float adc_samples[FFT_SIZE]; // ADC采样数据
// ADC中断服务程序
__interrupt void adc_isr(void)
{
static int sample_count = 0;
float sample;
// 读取ADC结果（假设12位ADC）
sample = (float)AdcRegs.ADCRESULT0 * 3.0 / 4095.0; // 转换为电压
// FIR滤波
sample = fir_filter(&fir_lpf, sample);
// IIR滤波
sample = iir_biquad_filter(&iir_biquad, sample);
// 添加到FFT缓冲区
realtime_fft_add_sample(&rt_fft, sample);
// 保存样本用于批处理
adc_samples[sample_count] = sample;
sample_count++;
if(sample_count >= FFT_SIZE)
{
sample_count = 0;
// 可以触发批处理FFT
}
// 清除ADC中断标志
AdcRegs.ADCST.bit.INT_SEQ1_CLR = 1;
PieCtrlRegs.PIEACK.all = PIEACK_GROUP1;
}
// 系统初始化
void System_Init(void)
{
// 初始化系统控制
InitSysCtrl();
// 初始化PIE控制
DINT;
InitPieCtrl();
IER = 0x0000;
IFR = 0x0000;
InitPieVectTable();
// 初始化ADC
InitAdc();
// 初始化FPU
InitFpu();
// 初始化ePWM（用于ADC触发）
InitEPwm();
}
// ADC配置
void ADC_Config(void)
{
// 配置ADC时钟
EALLOW;
SysCtrlRegs.HISPCP.all = ADC_MODCLK; // HSPCLK = SYSCLKOUT/ADC_MODCLK
EDIS;
// 配置ADC
AdcRegs.ADCTRL1.bit.ACQ_PS = 0x0F; // 采样窗口大小
AdcRegs.ADCTRL1.bit.CPS = 0; // 时钟预分频
AdcRegs.ADCTRL3.bit.ADCCLKPS = 0x3; // ADC内核时钟分频
AdcRegs.ADCTRL1.bit.SEQ_CASC = 1; // 级联序列模式
AdcRegs.ADCMAXCONV.bit.MAX_CONV = 0; // 1个转换
AdcRegs.ADCCHSELSEQ1.bit.CONV00 = 0; // 选择ADCINA0
// 使能中断
AdcRegs.ADCTRL2.bit.INT_ENA_SEQ1 = 1; // 使能SEQ1中断
AdcRegs.ADCTRL2.bit.INT_MOD_SEQ1 = 0; // 每个序列结束中断
// 配置中断
EALLOW;
PieVectTable.ADC = &adc_isr; // 设置中断向量
EDIS;
PieCtrlRegs.PIEIER1.bit.INTx6 = 1; // 使能ADC中断
IER |= M_INT1; // 使能CPU中断1
EINT; // 使能全局中断
}
// ePWM配置（用于ADC触发）
void EPWM_Config(void)
{
// 配置ePWM1
EPwm1Regs.TBPRD = 1500; // 周期 = 1500个TBCLK周期
EPwm1Regs.TBPHS.half.TBPHS = 0;
EPwm1Regs.TBCTL.bit.CTRMODE = TB_COUNT_UPDOWN; // 增减计数
// 配置ADC触发
EPwm1Regs.ETSEL.bit.SOCAEN = 1; // 使能SOCA
EPwm1Regs.ETSEL.bit.SOCASEL = 4; // 定时器等于0时触发
EPwm1Regs.ETPS.bit.SOCAPRD = 1; // 单次触发
// 启动ePWM
EPwm1Regs.TBCTL.bit.HSPCLKDIV = 0;
EPwm1Regs.TBCTL.bit.CLKDIV = 0;
}
int main(void)
{
// 系统初始化
System_Init();
// 初始化变量
int i;
// 生成窗函数
hanning_window(window, FFT_SIZE);
// 设计FIR滤波器（低通，截止频率100Hz，采样率1kHz）
design_fir_lpf(fir_coeffs, 65, 100.0, 1000.0);
// 初始化FIR滤波器
fir_init(&fir_lpf, fir_coeffs, fir_buffer, 65);
// 设计IIR滤波器（二阶低通，截止频率200Hz，采样率1kHz）
float fs = 1000.0;
float fc = 200.0;
float omega = 2.0 * PI * fc / fs;
float alpha = sin(omega) / (2.0 * 0.707); // Q=0.707
float b0 = (1.0 - cos(omega)) / 2.0;
float b1 = 1.0 - cos(omega);
float b2 = b0;
float a0 = 1.0 + alpha;
float a1 = -2.0 * cos(omega);
float a2 = 1.0 - alpha;
// 归一化
b0 /= a0;
b1 /= a0;
b2 /= a0;
a1 /= a0;
a2 /= a0;
// 初始化IIR滤波器
iir_biquad_init(&iir_biquad, b0, b1, b2, a1, a2);
// 初始化实时FFT
realtime_fft_init(&rt_fft);
// 配置ADC
ADC_Config();
// 配置ePWM
EPWM_Config();
// 主循环
while(1)
{
// 执行实时FFT处理
realtime_fft_process(&rt_fft);
// 查找主要频率分量
float max_mag = 0.0;
int max_index = 0;
float sampling_rate = 1000.0; // 1kHz采样率
for(i = 0; i < FFT_SIZE/2; i++)
{
if(rt_fft.mag[i] > max_mag)
{
max_mag = rt_fft.mag[i];
max_index = i;
}
}
// 计算频率
float freq = (float)max_index * sampling_rate / FFT_SIZE;
// 这里可以添加其他处理逻辑
// 例如：根据频率进行控制，或者通过通信接口发送结果
// 简单的延时
DELAY_US(1000);
}
return 0;
}
10. 批处理FFT和滤波示例
// 批处理示例
void batch_processing_example(void)
{
int i;
float filtered_data[FFT_SIZE];
float fft_real[FFT_SIZE];
float fft_imag[FFT_SIZE];
float magnitude[FFT_SIZE/2];
// 生成测试信号（50Hz正弦波 + 150Hz正弦波 + 噪声）
for(i = 0; i < FFT_SIZE; i++)
{
float t = (float)i / 1000.0; // 采样率1kHz
fft_input[i] = sin(2.0 * PI * 50.0 * t) +
0.5 * sin(2.0 * PI * 150.0 * t) +
0.1 * ((float)rand()/RAND_MAX - 0.5);
}
// FIR滤波
fir_filter_batch(&fir_lpf, fft_input, filtered_data, FFT_SIZE);
// 应用窗函数
apply_window(filtered_data, window, FFT_SIZE);
// 执行FFT
fft(filtered_data, fft_real, fft_imag, FFT_SIZE, LOG2_FFT_SIZE);
// 计算幅值谱
compute_magnitude(fft_real, fft_imag, magnitude, FFT_SIZE);
// 查找峰值频率
float max_mag = 0.0;
int max_index = 0;
float sampling_rate = 1000.0;
for(i = 1; i < FFT_SIZE/2; i++) // 跳过直流分量
{
if(magnitude[i] > max_mag)
{
max_mag = magnitude[i];
max_index = i;
}
}
float peak_freq = (float)max_index * sampling_rate / FFT_SIZE;
// 可以输出或使用peak_freq
}
11. 优化建议
// 使用DSP库进行优化
#include "DSP28x_FPU.h"
#include "FPU.h"
// 使用TI的FPU库进行复数乘法优化
#pragma CODE_SECTION(complex_multiply_fpu, "ramfuncs");
__interrupt void complex_multiply_fpu(Complex* a, Complex* b, Complex* result)
{
result->real = a->real * b->real - a->imag * b->imag;
result->imag = a->real * b->imag + a->imag * b->real;
}
// 使用查找表优化三角函数计算
float sin_table[360];
float cos_table[360];
void init_trig_tables(void)
{
int i;
for(i = 0; i < 360; i++)
{
float angle = (float)i * PI / 180.0;
sin_table[i] = sin(angle);
cos_table[i] = cos(angle);
}
}
// 优化蝶形运算
#pragma CODE_SECTION(butterfly_optimized, "ramfuncs");
void butterfly_optimized(Complex* data, int n, int m, float* sin_table, float* cos_table)
{
// 优化版本，使用查找表和循环展开
// ... 实现略 ...
}
参考代码 基于DSP 28335的傅里叶变化和滤波程序 www.youwenfan.com/contentcnv/71014.html
配置说明
- CMD文件配置：
MEMORY
{
PAGE 0: /* Program Memory */
PAGE 1: /* Data Memory */
/* 为FFT缓冲区分配特定内存区域 */
FFTBuffer : origin = 0x008000, length = 0x010000
FIRCoef : origin = 0x009000, length = 0x001000
IIRCoef : origin = 0x00A000, length = 0x001000
}
SECTIONS
{
/* 分配FFT缓冲区到特定区域 */
FFTBuffer : > FFTBuffer, PAGE = 1
FIRCoef : > FIRCoef, PAGE = 1
IIRCoef : > IIRCoef, PAGE = 1
}
-
编译选项：
- 启用FPU支持
- 优化级别设置为-O2或-O3
- 使用大内存模式
使用注意事项
- 实时性：FFT运算量较大，需要根据实际采样率和点数评估性能
- 内存使用：大点数FFT需要大量内存，合理分配内存区域
- 数值精度：定点实现时注意数据范围和精度
- 中断处理：FFT计算期间可能需要禁用中断以确保数据一致性
- 窗函数选择：根据应用需求选择合适的窗函数
