C语言进行离散傅里叶DFT变换~MATLAB验证
C语言进行离散傅里叶DFT变换~MATLAB验证
根据离散傅里叶变换的原始公式和自己编写复数计算函数进行离散傅里叶变换
对10000个点的加有噪声或干净的正弦波的数据进行离散傅里叶变换,生成10000个点的复数数据序列到文本文件中。
数据格式为实部+虚部,用空格或逗号隔开。
离散傅里叶变换的公式如下:
X
(
k
)
=
∑
n
=
0
N
?
1
x
(
n
)
exp
?
(
?
j
2
π
N
n
k
)
=
∑
n
=
0
N
?
1
x
(
n
)
W
N
n
k
\begin{aligned} X(k)&=\sum_{n=0}^{N-1}x(n)\exp(-j\frac{2\pi}{N}nk)\\ &=\sum_{n=0}^{N-1}x(n)W_N^nk \end{aligned}
X(k)?=n=0∑N?1?x(n)exp(?jN2π?nk)=n=0∑N?1?x(n)WNn?k?
带入欧拉公式
e
i
x
=
cos
?
x
+
i
(
sin
?
x
)
e^{ix}=\cos x+i(\sin x)
eix=cosx+i(sinx)
得到:
X
(
k
)
=
∑
n
=
0
N
?
1
x
(
n
)
cos
?
(
2
π
N
k
n
)
?
j
x
(
n
)
sin
?
(
2
π
N
k
n
)
X(k)=\sum_{n=0}^{N-1}x(n)\cos (\frac{2\pi}{N}kn)-jx(n)\sin (\frac{2\pi}{N}kn)
X(k)=n=0∑N?1?x(n)cos(N2π?kn)?jx(n)sin(N2π?kn)
幅值计算:
y
(
k
)
=
a
2
+
b
2
y(k)=\sqrt{a^2+b^2}
y(k)=a2+b2
?
#include
#include
#include
#define pi 3.14159
typedef struct
{
double real;//实部
double imag;//虚部
}complex;
//复数加法
void complex_add(complex a, complex b, complex *c)
{
c -> real = a.real + b.real;
c -> imag = a.imag + b.imag;
}
//复数乘法
void complex_mul(complex a, complex b, complex *c)
{
c -> real = a.real * b.real - a.imag * b.imag;
c -> imag = a.real * b.imag + a.imag * b.real;
}
//求模
double complex_mod(complex a)
{
double c = 0;
c = pow(a.real * a.real + a.imag + a.imag, 0.5);
return c;
}
int data_len(FILE*fp){
int len = 0,c;
while ((c = fgetc(fp)) != EOF)
if (c == '\n')
{
len++;
}
return len;
}
void Wn(int n1,int i,int k,complex *Wn) //定义旋转因子
{
//用欧拉公式分成实部和虚部
Wn->real=cos(2*pi*i*k/n1);
Wn->imag=-sin(2*pi*i*k/n1);
}
int main(int argc, char* argv[])
{
//检查命令行参数正确与否
if (argc != 4)
{
printf("用法:程序名 输入文件 输出文件 采样点数\n");
printf("程序名:编译后后的.exe名称\n");
printf("输入文件:所需进行计算的信号文件\n");
printf("输出文件:计算后输出的文件名\n");
printf("采样点数:所需计算的采样点数\n");
return -1;
}
FILE* fp1 = fopen(argv[1], "r");
int N = atoi(argv[2]);//定义采样点数
FILE* fp2 = fopen(argv[2], "w");
int len = 0;
double data = 0;
//检查文件能否被打开
if (fp1 == NULL)
{
printf("file error!\n");
return -1;
}
//检测文件行数
len = data_len(fp1);
complex x[len];//原序列
complex X[len];//DFT后的序列
double y[len];//幅值
printf("data len:%d\n",len);
rewind(fp1);
//拷贝文件内容
char* str = (char*)malloc(sizeof(char) * len);
for(int i = 0;i < len; i++)//把数据取出来
{
fscanf(fp1, "%s", str);
if (feof(fp1)) break;
data = atof(str);
x[i].real = data;
x[i].imag = 0;
}
for(int j =0;j
{//旋转因子计算
complex sum = {0,0};
for(int n = 0;n < len; n++)
{
complex wn,t;
Wn(len,n,j,&wn);
complex_mul(x[n],wn,&t);
complex_add(sum,t,&sum);
}
X[j].real = sum.real;
X[j].imag = sum.imag;
fprintf(fp2,"%lf %lf\n",X[j].real,X[j].imag);
}
// for(int k=0;k < len; k++)//幅值计算
// {
// double t = k/len;
// y[k] = sqrt(X[k].real * X[k].real
// + X[k].imag * X[k].imag);
// fprintf(fp2,"%lf %lf\n",t,y[k]);
// }
printf("success");
fclose(fp1);
fclose(fp2);
return 0;
}
在tcc中编译文件并运行
打开result1.txt,查看数据
将数据导入gnuplot中绘制频谱图
MATLAB中的代码如下
clear;
close all;
x = load('wave.txt');%读取数据
N = length(x);%检查矩阵的长度
a = zeros(N,1);%建立N行0矩阵,存放实部数据
b = zeros(N,1);%建立N行0矩阵,存放虚部数据
for k=0:N-1
for i=0:N-1
b(k+1)=b(k+1)+x(i+1)*sin((2*pi*k*i)/N);%虚部
a(k+1)=a(k+1)+x(i+1)*cos((2*pi*k*i)/N);%实部
end
c(k+1)=sqrt(a(k+1)^2+b(k+1)^2);%求幅度
end
t=(0:1:N-1);%时间数据
plot(t,c);%画频谱图
axis([0,300,0,5000]);
C语言进行离散傅里叶DFT变换~MATLAB验证相关教程
栈(简单介绍)
栈(简单介绍) 栈是一种只能在一端进行插入或者删除操作的线性表。表中允许插入,删除操作的一端叫做栈顶,另一端叫栈底。 但栈没有数据元素时称为空栈,插入数据操作为进栈或入栈,删除数据操作为出栈或退栈。 栈的特点为 后进先出 ,栈也成为后进先出表。如
2020-10-28 R语言统计与基础
2020-10-28 R语言统计与基础 Rstudio 比 R Gui 界面更友好,更适合新手使用 Rstudio界面介绍 Console 是一个命令行的窗口(运行的窗口),上方展开是 R 运行的脚本的窗口,在里面可以写很多命令 右上角是环境变量和历史的窗口 右下角,Plot 窗口是绘制图片的
VS 2019中利用C#语言在.Net Framework 4.5框架上开发简易串口发
VS 2019中利用C#语言在.Net Framework 4.5框架上开发简易串口发送数据软件 0. 前言 前置项目为:VS 2019中利用C#语言在.Net Framework 4.5框架上开发简易倒计时器 1. 项目实现 在Form1.cs[设计]中放入如下控件,并在对应的属性框中设置好属性。 拖入组件seria
15、HTML类
15、HTML类 对 HTML 进行分类(设置类),使我们能够为元素的类定义 CSS 样式。 为相同的类设置相同的样式,或者为不同的类设置不同的样式。 !DOCTYPE htmlhtmlheadstyle.cities { background-color:black; color:white; margin:20px; padding:20px;}/style/h
杭电 1096(C语言.while版)
杭电 1096(C语言.while版) 2020.10.28 我仅代表北京师范大学珠海分校海三C512的全体成员(当然不包括寿星)向 我们的舍友郑浩楠说一声生日快乐。 Input Input contains an integer N in the first line, and then N lines follow. Each line starts with a i
C语言编写程序:从键盘输入一个小写字母,该字母加密后变成其后
C语言编写程序:从键盘输入一个小写字母,该字母加密后变成其后继第2个字母输出。 C语言编写程序:从键盘输入一个小写字母,该字母加密后变成其后继第2个字母输出。 例如:a加密后变成c,b加密后变成d,z加密后变成b。 实验代码: #includestdio.hint main(){
进程
进程 进程: 进行就是正在进行中的程序,程序运行起来需要被加载到内存中。进程就是操作系统的描述,这个描述叫PCB(进程控制块),Linux下PCB有自己的名字叫task_struct。而操作系统就是使用task_struct结构体描述进程,使用双向链表来将这些结构体组织起来
C语言编写程序:从键盘输入三个数,将这三个数按从小到大排序并
C语言编写程序:从键盘输入三个数,将这三个数按从小到大排序并输出。 实验代码: #includestdio.hint main(){double x,y,z,t;printf(请从键盘上输入三个数:);scanf(%lf%lf%lf,x,y,z);if(xy){t=x;x=y;y=t;}if(xz){t=x;x=z;z=t;}if(yz){t=y;y=z;z=t;}printf(