📄 dftsave.c
字号:
/******************************************************
*
* file d:\cips\dftsave.c
*
* Functions: This file contains
* main
* dft
* invdft
* dft_2d
* invdft_2d
* odd_dft
* odd_invdft
* print_real_im
* print_2d_real_im
*
* Purpose:
* These functions perform forward and
* inverse discrete Fourier transforms
* in 1 and 2 dimensions. The basic algorithms
* are from "An Introduction to Digital
* Signal Processing," John H. Karl, Academic
* Press, 1989. The basic routines work for even
* dimensioned arrays. For odd dimensioned arrays
* you must use the odd routines.
*
* External Calls:
* none
*
* Modifications:
* 05 February 1991 - created
*
*****************************************************/
#include "d:\cips\cips.h"
#define N 5
#define M 5
#define pie 3.1425927
float box1[M][N] = {{0., 0., 0., 0., 0.},
{0.,-1.0 ,-1.0 ,-1.0 , 0.},
{0.,-1.0 , 9.0 ,-1.0 , 0.},
{0.,-1.0 ,-1.0 ,-1.0 , 0.},
{0., 0., 0., 0., 0.}};
float box2[M][N] = {{0., 0., 0., 0., 0.},
{0., 0., 0., 0., 0.},
{0., 0., 0., 0., 0.},
{0., 0., 0., 0., 0.},
{0., 0., 0., 0., 0.}};
float box1_out[M][N] = {{3.5, 3.0, 2.5, 3.0, 3.5},
{3.0, 2.5, 1.7, 2.5, 3.0},
{2.8, 1.7, 1.0, 1.7, 2.8},
{3.0, 2.5 , 1.7, 2.5, 3.0},
{3.5, 3.0, 2.5, 3.0, 3.5}};
float box2_out[M][N] = {{0., 0., 0., 0., 0.},
{0., 0., 0., 0., 0.},
{0., 0., 0., 0., 0.},
{0., 0., 0., 0., 0.},
{0., 0., 0., 0., 0.}};
float x[N] = {1., 1., 1.};
float y[N] = {0., 0., 0.};
float r[N], i[N];
main()
{
/*** 2-D inverse DFT ***/
printf("\nDFT> This is the 2D input in freq domain\n");
print_2d_real_im(box1_out, box2_out);
printf("\nDFT> Now take the inverse DFT and look at time domain\n");
odd_invdft_2d(box1_out, box2_out, box1, box2);
print_2d_real_im(box1, box2);
/**** 1-D forward then inverse DFT
printf("\nDFT> This is the input\n");
print_real_im(x, y);
dft(x, y, r, i);
printf("\nDFT> This is the output\n");
print_real_im(r, i);
invdft(r, i, x, y);
printf("\nDFT> Now reverse the process and look at the data\n");
print_real_im(x, y);
****/
/*** 2-D forward and inverse DFT
printf("\nDFT> This is the 2D input\n");
print_2d_real_im(box1, box2);
dft_2d(box1, box2, box1_out, box2_out);
printf("\nDFT> This is the 2D output\n");
print_2d_real_im(box1_out, box2_out);
printf("\nDFT> Now reverse the process and look at the data\n");
invdft_2d(box1_out, box2_out, box1, box2);
print_2d_real_im(box1, box2);
****/
/*** 2-D inverse DFT
printf("\nDFT> This is the 2D input in freq domain\n");
print_2d_real_im(box1_out, box2_out);
printf("\nDFT> Now take the inverse DFT and look at time domain\n");
invdft_2d(box1_out, box2_out, box1, box2);
print_2d_real_im(box1, box2);
*****/
} /* ends main */
/*
This is a simple print routine to
look at the 2D real and imaginary
for small M N
*/
print_2d_real_im(a, b)
float a[M][N], b[M][N];
{
int i, k;
printf("\nDFT> ");
for(i=0; i<M; i++){
printf("\nDFT> ");
for(k=0; k<N; k++){
printf(" %2.4f+j%2.4f", a[i][k], b[i][k]);
}
}
} /* ends print_2d_real_im */
/*
This is a simple print routine to
look at the 1D real and imaginary
for small N
*/
print_real_im(a, b)
float a[], b[];
{
int i;
printf("\nDFT> ");
for(i=0; i<N; i++)
printf("\nDFT> %f + j%f", a[i], b[i]);
} /* ends print_real_im */
/*
This is the 1D forward DFT.
This is the centered format.
This runs from -N/2 to N/2 + 1.
This works for N even.
*/
dft(x, y, r, i)
float x[], y[], r[], i[];
{
int k, j;
float c, s, p, w, x0, y0;
for(k=-N/2; k<=N/2 -1; k++){
w = 2. * pie * k/N;
x0 = 0;
y0 = 0;
for(j=-N/2; j<=N/2 - 1; j++){
/*p = w * j;*/
p = 2. * pie * k * j/N;
c = cos(p);
s = sin(p);
x0 = x0 + c*x[j+N/2] + s*y[j+N/2];
y0 = y0 + c*y[j+N/2] - s*x[j+N/2];
} /* ends loop over j */
r[k+N/2] = x0;
i[k+N/2] = y0;
} /* ends loop over k */
} /* ends dft */
/*
This is the 1D reverse DFT.
This is the centered format.
This runs from -N/2 to N/2 + 1.
This works for N even.
*/
invdft(x, y, r, i)
float x[], y[], r[], i[];
{
int k, j;
float c, s, p, w, x0, y0;
for(k=-N/2; k<=N/2 -1; k++){
w = -1. * 2. * pie * k/N;
x0 = 0;
y0 = 0;
for(j=-N/2; j<=N/2 - 1; j++){
p = w * j;
c = cos(p);
s = sin(p);
x0 = x0 + c*x[j+N/2] + s*y[j+N/2];
y0 = y0 + c*y[j+N/2] - s*x[j+N/2];
} /* ends loop over j */
r[k+N/2] = x0/N;
i[k+N/2] = y0/N;
} /* ends loop over k */
} /* ends invdft */
/*
This is the 1D forward DFT.
This is the centered format.
This runs from -N/2 to N/2.
This works for N odd.
*/
odd_dft(x, y, r, i)
float x[], y[], r[], i[];
{
int k, j;
float c, s, p, w, x0, y0;
for(k=-N/2; k<=N/2; k++){
w = 2. * pie * k/N;
x0 = 0;
y0 = 0;
for(j=-N/2; j<=N/2; j++){
/*p = w * j;*/
p = 2. * pie * k * j/N;
c = cos(p);
s = sin(p);
x0 = x0 + c*x[j+N/2] + s*y[j+N/2];
y0 = y0 + c*y[j+N/2] - s*x[j+N/2];
} /* ends loop over j */
r[k+N/2] = x0;
i[k+N/2] = y0;
} /* ends loop over k */
} /* ends odd_dft */
/*
This is the 1D reverse DFT.
This is the centered format.
This runs from -N/2 to N/2.
This works for N odd.
*/
odd_invdft(x, y, r, i)
float x[], y[], r[], i[];
{
int k, j;
float c, s, p, w, x0, y0;
for(k=-N/2; k<=N/2; k++){
w = -1. * 2. * pie * k/N;
x0 = 0;
y0 = 0;
for(j=-N/2; j<=N/2; j++){
p = w * j;
c = cos(p);
s = sin(p);
x0 = x0 + c*x[j+N/2] + s*y[j+N/2];
y0 = y0 + c*y[j+N/2] - s*x[j+N/2];
} /* ends loop over j */
r[k+N/2] = x0/N;
i[k+N/2] = y0/N;
} /* ends loop over k */
} /* ends odd_invdft */
/*
This is the forward 2D DFT.
This is the centered format.
This runs from -N/2 to N/2 + 1
This runs from -M/2 to M/2 + 1
This works for MxN matrices
where M and N are even.
x is the input array real part.
y is the input array imaginary part.
r is the output array real part.
i is the output array imaginary part.
The time domain h has subscripts
m n. The freq domain H has
subscripts u v. These are both
varied over M N.
The angle p = -2jpienu/N - 2jpiemv/M
p = 2jpie[ (-nuM) - (mvN)]/MN
*/
dft_2d(x, y, r, i)
float x[M][N], y[M][N], r[M][N], i[M][N];
{
int n, m, u, v;
float c, s, p, w, x0, y0;
for(v=-M/2; v<=M/2 -1; v++){
for(u=-N/2; u<=N/2 -1; u++){
x0 = 0;
y0 = 0;
for(m=-M/2; m<=M/2 - 1; m++){
for(n=-N/2; n<=N/2 - 1; n++){
/* you can probably separate the following
to increase speed */
p = 2. * pie * (n*u*M + m*v*N) / (N*M);
c = cos(p);
s = sin(p);
x0 = x0 + c*x[m+M/2][n+N/2] + s*y[m+M/2][n+N/2];
y0 = y0 + c*y[m+M/2][n+N/2] - s*x[m+M/2][n+N/2];
} /* ends loop over n */
} /* ends loop over m */
r[v+M/2][u+N/2] = x0;
i[v+M/2][u+N/2] = y0;
} /* ends loop over u */
} /* ends loop over v */
} /* ends dft_2d */
/*
This is the reverse 2D DFT.
This is the centered format.
This runs from -N/2 to N/2 + 1
This runs from -M/2 to M/2 + 1
This works for MxN matrices
where M and N are even.
x is the input array real part.
y is the input array imaginary part.
r is the output array real part.
i is the output array imaginary part.
The time domain h has subscripts
m n. The freq domain H has
subscripts u v. These are both
varied over M N.
The angle p = -2jpienu/N - 2jpiemv/M
p = 2jpie[ (-nuM) - (mvN)]/MN
*/
invdft_2d(x, y, r, i)
float x[M][N], y[M][N], r[M][N], i[M][N];
{
int n, m, u, v;
float c, s, p, w, x0, y0;
for(v=-M/2; v<=M/2 -1; v++){
for(u=-N/2; u<=N/2 -1; u++){
x0 = 0;
y0 = 0;
for(m=-M/2; m<=M/2 - 1; m++){
for(n=-N/2; n<=N/2 - 1; n++){
/* you can probably separate the following
to increase speed */
p = 2. * pie * (-1*n*u*M - m*v*N) / (N*M);
c = cos(p);
s = sin(p);
x0 = x0 + c*x[m+M/2][n+N/2] + s*y[m+M/2][n+N/2];
y0 = y0 + c*y[m+M/2][n+N/2] - s*x[m+M/2][n+N/2];
} /* ends loop over n */
} /* ends loop over m */
r[v+M/2][u+N/2] = x0/(M*N);
i[v+M/2][u+N/2] = y0/(M*N);
} /* ends loop over u */
} /* ends loop over v */
} /* ends invdft_2d */
/*
This is the forward 2D DFT.
This is the centered format.
This runs from -N/2 to N/2
This runs from -M/2 to M/2
This works for MxN matrices
where M and N are odd.
x is the input array real part.
y is the input array imaginary part.
r is the output array real part.
i is the output array imaginary part.
The time domain h has subscripts
m n. The freq domain H has
subscripts u v. These are both
varied over M N.
The angle p = -2jpienu/N - 2jpiemv/M
p = 2jpie[ (-nuM) - (mvN)]/MN
*/
odd_dft_2d(x, y, r, i)
float x[M][N], y[M][N], r[M][N], i[M][N];
{
int n, m, u, v;
float c, s, p, w, x0, y0;
for(v=-M/2; v<=M/2; v++){
for(u=-N/2; u<=N/2; u++){
x0 = 0;
y0 = 0;
for(m=-M/2; m<=M/2; m++){
for(n=-N/2; n<=N/2; n++){
/* you can probably separate the following
to increase speed */
p = 2. * pie * (n*u*M + m*v*N) / (N*M);
c = cos(p);
s = sin(p);
x0 = x0 + c*x[m+M/2][n+N/2] + s*y[m+M/2][n+N/2];
y0 = y0 + c*y[m+M/2][n+N/2] - s*x[m+M/2][n+N/2];
} /* ends loop over n */
} /* ends loop over m */
r[v+M/2][u+N/2] = x0;
i[v+M/2][u+N/2] = y0;
} /* ends loop over u */
} /* ends loop over v */
} /* ends odd_dft_2d */
/*
This is the reverse 2D DFT.
This is the centered format.
This runs from -N/2 to N/2
This runs from -M/2 to M/2
This works for MxN matrices
where M and N are odd.
x is the input array real part.
y is the input array imaginary part.
r is the output array real part.
i is the output array imaginary part.
The time domain h has subscripts
m n. The freq domain H has
subscripts u v. These are both
varied over M N.
The angle p = -2jpienu/N - 2jpiemv/M
p = 2jpie[ (-nuM) - (mvN)]/MN
*/
odd_invdft_2d(x, y, r, i)
float x[M][N], y[M][N], r[M][N], i[M][N];
{
int n, m, u, v;
float c, s, p, w, x0, y0;
for(v=-M/2; v<=M/2; v++){
for(u=-N/2; u<=N/2; u++){
x0 = 0;
y0 = 0;
for(m=-M/2; m<=M/2; m++){
for(n=-N/2; n<=N/2; n++){
/* you can probably separate the following
to increase speed */
p = 2. * pie * (-1*n*u*M - m*v*N) / (N*M);
c = cos(p);
s = sin(p);
x0 = x0 + c*x[m+M/2][n+N/2] + s*y[m+M/2][n+N/2];
y0 = y0 + c*y[m+M/2][n+N/2] - s*x[m+M/2][n+N/2];
} /* ends loop over n */
} /* ends loop over m */
r[v+M/2][u+N/2] = x0/(M*N);
i[v+M/2][u+N/2] = y0/(M*N);
} /* ends loop over u */
} /* ends loop over v */
} /* ends odd_invdft_2d */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -