我正在尝试为2D FFT编写C程序。我们的想法是传递具有预分配内存的结构指针,其中函数将修改这些内存地址。
complex_num是一种支持复数运算的结构,如加法,减法,乘法,除法,幂,复指数。我放置printf的地方(“Hello \ n”);语句是我的代码中发生分段错误的地方。谢谢你的期待。
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
typedef struct complex_num_struct {
double real;
double imag;
} complex_num;
complex_num add(complex_num num1, complex_num num2);
complex_num subtract(complex_num num1, complex_num num2);
complex_num multiply(complex_num num1, complex_num num2);
complex_num divide(complex_num num1, complex_num num2);
complex_num power(complex_num num, double n);
complex_num complex_exp(double theta);
double magnitude(complex_num num);
void display_complex_matrix(complex_num** matrix, int height, int width);
void fft2(complex_num** input, complex_num** output, int height, int width);
void fft_driver(complex_num input[], complex_num output[], int n, int step);
void fft(complex_num input[], complex_num output[], int n);
void display_complex_vec(complex_num * vec, int n);
complex_num conjugate(complex_num num);
complex_num ** matrix;
complex_num ** temp;
void main() {
complex_num one;
one.real = 1;
one.imag = 0;
complex_num height, width;
height.real = 4;
width.real = 4;
height.imag = 0;
width.imag = 0;
complex_num zero;
zero.real = 0;
zero.imag = 0;
complex_num input1[] = {one, one, one, one};
complex_num input2[] = {one, one, one, one};
complex_num input3[] = {zero, zero, zero, zero};
complex_num input4[] = {zero, zero, zero, zero};
complex_num* input[] = {input1, input2, input3, input4};
complex_num* output[] = {input1, input2, input3, input4};
complex_num* conj[] = {input1, input2, input3, input4};
fft2(input, output, 4, 4);
display_complex_matrix(output, 4, 4);
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
conj[i][j] = conjugate(output[i][j]);
}
}
fft2(conj, output, 4, 4);
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
output[i][j] = conjugate(divide(conj[i][j], (multiply(height, width))));
}
}
display_complex_matrix(output, 4, 4);
}
void fft2(complex_num** input, complex_num** output, int height, int width) {
matrix = (complex_num **) malloc(height * sizeof(complex_num*));
temp = (complex_num **) malloc(width * sizeof(complex_num*));
for (int i = 0; i < height; ++i) {
fft(input[i], temp[i], width);
}
for (int i = 0; i < height; ++i) {
//matrix[i] = (complex_num*) malloc(height * sizeof(complex_num));
for (int j = 0; j < width; ++j) {
//matrix[i][j] = temp[j][i];
printf("%lf", temp[j][i].real);
}
}
for (int i = 0; i < height; ++i) {
fft(matrix[i], temp[i], width);
}
for (int i = 0; i < height; ++i) {
for (int j = 0; j < width; ++j) {
output[i][j] = temp[j][i];
}
}
}
void display_complex_vec(complex_num vec[], int n) {
printf("[");
for (int i = 0; i < n; ++i) {
printf("(%g, %g), ", vec[i].real, vec[i].imag);
}
printf("]\n");
}
void display_complex_matrix(complex_num** matrix, int height, int width) {
for (int i = 0; i < height; ++i) {
printf("[");
for (int j = 0; j < width; ++j) {
printf("(%g, %g), ", matrix[i][j].real, matrix[i][j].imag);
}
printf("]\n");
}
}
void fft_driver(complex_num input[], complex_num output[], int n, int step) {
complex_num diff, cexp;
if (step < n) {
fft_driver(output, input, n, step * 2);
fft_driver(output + step, input + step, n, step * 2);
for (int i = 0; i < n; i += 2 * step) {
cexp = complex_exp(-M_PI * (double) i / (double) n);
diff = multiply(complex_exp(-M_PI * (double) i / (double) n), output[i + step]);
input[i / 2] = add(output[i], diff);
input[(i + n) / 2] = subtract(output[i], diff);
//printf("(%g, %g) * (%g, %g)\n", input[i / 2].real, input[i / 2].imag, input[(i + n) / 2].real, input[(i + n) / 2].imag);
}
}
}
void fft(complex_num* input, complex_num* output, int n) {
complex_num* dummy = (complex_num *)malloc(n * sizeof(complex_num));
complex_num* inp = (complex_num *)malloc(n * sizeof(complex_num));
for (int i = 0; i < n; ++i) {
dummy[i].real = 0;
dummy[i].imag = 0;
inp[i].real = 0;
inp[i].imag = 0;
}
printf("%lf", inp[1].real);
for (int i = 0; i < n; ++i) {
dummy[i] = input[i];
inp[i] = input[i];
}
fft_driver(inp, dummy, n, 1);
printf("Hello\n");
for (int i = 0; i < n; ++i) {
output[i] = inp[i];
}
}
complex_num complex_exp(double theta) {
complex_num result;
result.real = cos(theta);
result.imag = sin(theta);
return result;
}
double magnitude(complex_num num) {
return sqrt(pow(num.real, 2) + pow(num.imag, 2));
}
complex_num add(complex_num num1, complex_num num2) {
complex_num result;
result.real = num1.real + num2.real;
result.imag = num1.imag + num2.imag;
return result;
}
complex_num conjugate(complex_num num) {
complex_num result;
result.real = num.real;
result.imag = -num.imag;
return result;
}
complex_num subtract(complex_num num1, complex_num num2) {
complex_num result;
result.real = num1.real - num2.real;
result.imag = num1.imag - num2.imag;
return result;
}
complex_num multiply(complex_num num1, complex_num num2) {
complex_num result;
result.real = num1.real * num2.real - num1.imag * num2.imag;
result.imag = num1.real * num2.imag + num1.imag * num2.real;
return result;
}
complex_num divide(complex_num num1, complex_num num2) {
complex_num result;
double magnitude_square = pow(num2.real, 2) + pow(num2.imag, 2);
num2.imag = -num2.imag;
result = multiply(num1, num2);
result.real = result.real / magnitude_square;
result.imag = result.imag / magnitude_square;
return result;
}
complex_num power(complex_num num, double n) {
complex_num result;
double magnitude = sqrt(pow(num.real, 2) + pow(num.imag, 2));
result.real = pow(magnitude, n) * cos(n * atan(num.imag / num.real));
result.imag = pow(magnitude, n) * sin(n * atan(num.imag / num.real));
return result;
}
上面的代码产生以下输出:
0.000000Hello
Segmentation fault (core dumped)
在调试时,我发现分段错误发生在printf语句之后的循环中我的代码的以下函数中:
void fft(complex_num* input, complex_num* output, int n) {
complex_num* dummy = (complex_num *)malloc(n * sizeof(complex_num));
complex_num* inp = (complex_num *)malloc(n * sizeof(complex_num));
for (int i = 0; i < n; ++i) {
dummy[i].real = 0;
dummy[i].imag = 0;
inp[i].real = 0;
inp[i].imag = 0;
}
printf("%lf", inp[1].real);
for (int i = 0; i < n; ++i) {
dummy[i] = input[i];
inp[i] = input[i];
}
fft_driver(inp, dummy, n, 1);
printf("Hello\n");
for (int i = 0; i < n; ++i) {
output[i] = inp[i];
}
}
所需的输出应如下: -
0.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
4.0000004.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
[(8, 0), (0, 0), (0, 0), (0, 0), ]
[(4, -4), (0, 0), (0, 0), (0, 0), ]
[(0, 0), (0, 0), (0, 0), (0, 0), ]
[(4, 4), (0, 0), (0, 0), (0, 0), ]
0.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
4.0000004.0000004.0000004.0000004.0000004.0000004.0000004.0000004.0000004.0000004.0000004.0000004.0000004.0000004.0000004.0000000.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
[(1, 1), (1, 1), (1, 1), (1, 1), ]
[(0.0625, -0.0625), (0.0625, -0.0625), (0.0625, -0.0625), (0.0625, -0.0625), ]
[(0.0625, -0.0625), (0.0625, -0.0625), (0.0625, -0.0625), (0.0625, -0.0625), ]
[(0.0625, -0.0625), (0.0625, -0.0625), (0.0625, -0.0625), (0.0625, -0.0625), ]
您有一个分段错误,因为您在void fft(complex_num* input, complex_num* output, int n)
的一行未写入的地址写入:
output[i] = inp[i];
输出来自行中的fft2
fft(input[i], temp[i], width);
temp设定的地方
temp = (complex_num **) malloc(width * sizeof(complex_num*));
但是temp[i]
没有初始化
也许fft(input[i], temp[i], width);
应该是fft(input[i], &temp[i], width);
而fft成为:
void fft(complex_num* input, complex_num** output, int n) {
complex_num* dummy = (complex_num *)malloc(n * sizeof(complex_num));
complex_num* inp = (complex_num *)malloc(n * sizeof(complex_num));
for (int i = 0; i < n; ++i) {
dummy[i].real = 0;
dummy[i].imag = 0;
inp[i].real = 0;
inp[i].imag = 0;
}
printf("%lf", inp[1].real);
for (int i = 0; i < n; ++i) {
dummy[i] = input[i];
inp[i] = input[i];
}
fft_driver(inp, dummy, n, 1);
printf("Hello\n");
*output = inp;
}
无论如何,在fft2中还有另一个问题:
fft(matrix[i], &temp[i], width);
因为matrix[i]
的初始化被注释,matrix[i]
没有被初始化,并且在fft中你正在访问这些行中的未知地址
dummy[i] = input[i];
inp[i] = input[i];
行//matrix[i] = (complex_num*) malloc(height * sizeof(complex_num));
不得在评论中,但也要将分配的大小更改为:
matrix[i] = (complex_num*) malloc(width * sizeof(complex_num));
在matrix[i][j] = temp[j][i];
以下的行不得发表评论
执行所有更改(我使用3.1415927表示M_PI)
/tmp % ./a.out
0.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
4.0000004.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
[(8, 0), (0, 0), (0, 0), (0, 0), ]
[(4, -4), (0, 0), (0, 0), (0, 0), ]
[(0, 0), (0, 0), (0, 0), (0, 0), ]
[(4, 4), (0, 0), (0, 0), (0, 0), ]
0.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
8.0000004.0000000.0000004.0000008.0000004.0000000.0000004.0000008.0000004.0000000.0000004.0000008.0000004.0000000.0000004.0000000.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
[(1, -0), (1, -0), (1, -0), (1, -0), ]
[(1, -2.8606e-18), (1, -2.8606e-18), (1, -2.8606e-18), (1, -2.8606e-18), ]
[(0, 0), (0, -0), (0, -0), (0, -0), ]
[(0, 2.8606e-18), (0, 2.8606e-18), (0, 2.8606e-18), (0, 2.8606e-18), ]
我不知道这是否是预期的输出,但至少valgrind没有检测到任何非法的内存访问,也没有检测到非初始化的值。
你有内存泄漏,删除它们:
在fft2替换
for (int i = 0; i < height; ++i) {
fft(matrix[i], &temp[i], width);
}
通过
for (int i = 0; i < height; ++i) {
free(temp[i]);
fft(matrix[i], &temp[i], width);
free(matrix[i]);
}
free(matrix);
并在最后添加
for (int i = 0; i < height; ++i)
free(temp[i]);
free(temp);
在fft结束前添加free(dummy);
这些更改会删除所有内存泄漏:
/tmp % valgrind --leak-check=full ./a.out
==12924== Memcheck, a memory error detector
==12924== Copyright (C) 2002-2012, and GNU GPL'd, by Julian Seward et al.
==12924== Using Valgrind-3.8.1 and LibVEX; rerun with -h for copyright info
==12924== Command: ./a.out
==12924==
0.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
4.0000004.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.0000000.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
[(8, 0), (0, 0), (0, 0), (0, 0), ]
[(4, -4), (0, 0), (0, 0), (0, 0), ]
[(0, 0), (0, 0), (0, 0), (0, 0), ]
[(4, 4), (0, 0), (0, 0), (0, 0), ]
0.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
8.0000004.0000000.0000004.0000008.0000004.0000000.0000004.0000008.0000004.0000000.0000004.0000008.0000004.0000000.0000004.0000000.000000Hello
0.000000Hello
0.000000Hello
0.000000Hello
[(1, -0), (1, -0), (1, -0), (1, -0), ]
[(1, -2.8606e-18), (1, -2.8606e-18), (1, -2.8606e-18), (1, -2.8606e-18), ]
[(0, 0), (0, -0), (0, -0), (0, -0), ]
[(0, 2.8606e-18), (0, 2.8606e-18), (0, 2.8606e-18), (0, 2.8606e-18), ]
==12924==
==12924== HEAP SUMMARY:
==12924== in use at exit: 0 bytes in 0 blocks
==12924== total heap usage: 44 allocs, 44 frees, 2,688 bytes allocated
==12924==
==12924== All heap blocks were freed -- no leaks are possible
==12924==
==12924== For counts of detected and suppressed errors, rerun with: -v
==12924== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 6 from 6)
注意主要必须是int main()
而不是void main()
如果我把所有代码都用来帮助你,因为有很多变化:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#ifndef M_PI
#define M_PI 3.1415927
#endif
typedef struct complex_num_struct {
double real;
double imag;
} complex_num;
complex_num add(complex_num num1, complex_num num2);
complex_num subtract(complex_num num1, complex_num num2);
complex_num multiply(complex_num num1, complex_num num2);
complex_num divide(complex_num num1, complex_num num2);
complex_num power(complex_num num, double n);
complex_num complex_exp(double theta);
double magnitude(complex_num num);
void display_complex_matrix(complex_num** matrix, int height, int width);
void fft2(complex_num** input, complex_num** output, int height, int width);
void fft_driver(complex_num input[], complex_num output[], int n, int step);
void fft(complex_num input[], complex_num * output[], int n);
void display_complex_vec(complex_num * vec, int n);
complex_num conjugate(complex_num num);
complex_num ** matrix;
complex_num ** temp;
int main() {
complex_num one;
one.real = 1;
one.imag = 0;
complex_num height, width;
height.real = 4;
width.real = 4;
height.imag = 0;
width.imag = 0;
complex_num zero;
zero.real = 0;
zero.imag = 0;
complex_num input1[] = {one, one, one, one};
complex_num input2[] = {one, one, one, one};
complex_num input3[] = {zero, zero, zero, zero};
complex_num input4[] = {zero, zero, zero, zero};
complex_num* input[] = {input1, input2, input3, input4};
complex_num* output[] = {input1, input2, input3, input4};
complex_num* conj[] = {input1, input2, input3, input4};
fft2(input, output, 4, 4);
display_complex_matrix(output, 4, 4);
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
conj[i][j] = conjugate(output[i][j]);
}
}
fft2(conj, output, 4, 4);
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 4; ++j) {
output[i][j] = conjugate(divide(conj[i][j], (multiply(height, width))));
}
}
display_complex_matrix(output, 4, 4);
}
void fft2(complex_num** input, complex_num** output, int height, int width) {
matrix = (complex_num **) malloc(height * sizeof(complex_num*));
temp = (complex_num **) malloc(width * sizeof(complex_num*));
for (int i = 0; i < height; ++i) {
fft(input[i], &temp[i], width);
}
for (int i = 0; i < height; ++i) {
matrix[i] = (complex_num*) malloc(width * sizeof(complex_num));
for (int j = 0; j < width; ++j) {
matrix[i][j] = temp[j][i];
printf("%lf", temp[j][i].real);
}
}
for (int i = 0; i < height; ++i) {
free(temp[i]);
fft(matrix[i], &temp[i], width);
free(matrix[i]);
}
free(matrix);
for (int i = 0; i < height; ++i) {
for (int j = 0; j < width; ++j) {
output[i][j] = temp[j][i];
}
}
for (int i = 0; i < height; ++i)
free(temp[i]);
free(temp);
}
void display_complex_vec(complex_num vec[], int n) {
printf("[");
for (int i = 0; i < n; ++i) {
printf("(%g, %g), ", vec[i].real, vec[i].imag);
}
printf("]\n");
}
void display_complex_matrix(complex_num** matrix, int height, int width) {
for (int i = 0; i < height; ++i) {
printf("[");
for (int j = 0; j < width; ++j) {
printf("(%g, %g), ", matrix[i][j].real, matrix[i][j].imag);
}
printf("]\n");
}
}
void fft_driver(complex_num input[], complex_num output[], int n, int step) {
complex_num diff, cexp;
if (step < n) {
fft_driver(output, input, n, step * 2);
fft_driver(output + step, input + step, n, step * 2);
for (int i = 0; i < n; i += 2 * step) {
cexp = complex_exp(-M_PI * (double) i / (double) n);
diff = multiply(complex_exp(-M_PI * (double) i / (double) n), output[i + step]);
input[i / 2] = add(output[i], diff);
input[(i + n) / 2] = subtract(output[i], diff);
//printf("(%g, %g) * (%g, %g)\n", input[i / 2].real, input[i / 2].imag, input[(i + n) / 2].real, input[(i + n) / 2].imag);
}
}
}
void fft(complex_num* input, complex_num** output, int n) {
complex_num* dummy = (complex_num *)malloc(n * sizeof(complex_num));
complex_num* inp = (complex_num *)malloc(n * sizeof(complex_num));
for (int i = 0; i < n; ++i) {
dummy[i].real = 0;
dummy[i].imag = 0;
inp[i].real = 0;
inp[i].imag = 0;
}
printf("%lf", inp[1].real);
for (int i = 0; i < n; ++i) {
dummy[i] = input[i];
inp[i] = input[i];
}
fft_driver(inp, dummy, n, 1);
free(dummy);
printf("Hello\n");
*output = inp;
}
complex_num complex_exp(double theta) {
complex_num result;
result.real = cos(theta);
result.imag = sin(theta);
return result;
}
double magnitude(complex_num num) {
return sqrt(pow(num.real, 2) + pow(num.imag, 2));
}
complex_num add(complex_num num1, complex_num num2) {
complex_num result;
result.real = num1.real + num2.real;
result.imag = num1.imag + num2.imag;
return result;
}
complex_num conjugate(complex_num num) {
complex_num result;
result.real = num.real;
result.imag = -num.imag;
return result;
}
complex_num subtract(complex_num num1, complex_num num2) {
complex_num result;
result.real = num1.real - num2.real;
result.imag = num1.imag - num2.imag;
return result;
}
complex_num multiply(complex_num num1, complex_num num2) {
complex_num result;
result.real = num1.real * num2.real - num1.imag * num2.imag;
result.imag = num1.real * num2.imag + num1.imag * num2.real;
return result;
}
complex_num divide(complex_num num1, complex_num num2) {
complex_num result;
double magnitude_square = pow(num2.real, 2) + pow(num2.imag, 2);
num2.imag = -num2.imag;
result = multiply(num1, num2);
result.real = result.real / magnitude_square;
result.imag = result.imag / magnitude_square;
return result;
}
complex_num power(complex_num num, double n) {
complex_num result;
double magnitude = sqrt(pow(num.real, 2) + pow(num.imag, 2));
result.real = pow(magnitude, n) * cos(n * atan(num.imag / num.real));
result.imag = pow(magnitude, n) * sin(n * atan(num.imag / num.real));
return result;
}