信号与系统实验2：信号的矩形脉冲抽样与恢复

这是信号与系统的实验报告。我信号与系统总共就扣了几分，也不知道是期末扣的，还是实验扣的。即使全是实验扣的，这实验也算做的挺好的。于是把报告发出来给大家参考。 这是第二个实验，信号的矩形脉冲抽样与恢复

一、实验结果展示

各频域图像

  1. \(F(\omega)\)原始频域图像

    image-20220604130549021

  2. \(0.2\text{Hz}\)采样频域图像

    image-20220604130602889

  3. \(0.5\text{Hz}\)采样频域图像

    image-20220604130625486

  4. \(1\text{Hz}\)采样频域图像：

    image-20220604130636810

  5. \(0.2\text{Hz}\)采样滤波后频域图像：

    image-20220604130648381

  6. \(0.5\text{Hz}\)采样滤波后频域图像：

    image-20220604130659596

  7. \(1\text{Hz}\)采样滤波后频域图像：

    image-20220604130713220

各时域图像

  1. \(F(\omega)\)傅里叶反变换后时域图像：

    image-20220604130724746

  2. \(0.2\text{Hz}\)采样滤波后恢复的时域图像

    image-20220604130735531

  3. \(0.5\text{Hz}\)采样滤波后恢复的时域图像：

    image-20220604130746852

  4. \(1\text{Hz}\)采样滤波后恢复的时域图像

    image-20220604130758457

二、实验结果分析

由奈奎斯特抽样定理，一个带限信号\(f(t)\),如果其频谱存在在频域（角频率）区间\([-\omega_m,\omega_m]\)，则可用抽样值唯一表示\(f(t)\),抽样值的间隔不能大于\(T_s=\frac{1}{2f_m}\),其中\(f_m=\frac{\omega_m}{2\pi}\).

在本实验中，\(\omega_m=\frac \pi2\)，即采样的频率需要达到\(2\frac{\omega_m}{2\pi}=0.5\text{Hz}\)。所以\(0.5\text{Hz}\)和\(1\text{Hz}\)采样后可以恢复，而\(0.2\text{Hz}\)采样后就不能恢复出原来的波形。

三、实验代码

  1. 实验用代码

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    #include "arrayToPPM.h"
    #define maxn 5005

    const double eps = 1e-11;
    const double PI = 3.1415926535;
    typedef double (*fun_p)(double);

    double F(double omega) {
    	return ((omega >= -0.5 * PI && omega <= 0.5 * PI) ? cos(omega) : 0);
    }

    double fourierTransform(double w, fun_p f, double left, double right, double step) {
    	double ans = 0.0;
    	for (double t = left; t <= right; t += step) {
    		ans = ans + f(t) * cos(w * t) * step;
    	}
    	return ans;
    }

    double fourierInvTransform(double t, fun_p f, double left, double right, double step) {
    	double ans = 0.0;
    	for (double w = left; w <= right; w += step) {
    		ans = ans + f(w) * cos(w * t) * step;
    	}
    	return ans / (2.0 * PI);
    }

    double sample(double t, double fre, double tau) { //时域矩形抽样函数
    	double clc = 1.0 / fre;
    	t = fabs(t);
    	double res = fmod(t, clc);
    	if (res > (clc / 2.0)) res -= clc;
    	return (fabs(res) <= (tau / 2.0) ? 1.0 : 0.0);
    }

    double filter(double w, double fre, double tau) {
    	if (w >= -PI / 2.0 && w <= PI / 2.0)
    		return (1.0 / fre) / tau;
    	else
    		return 0.0;
    }


    int main() {
    	double t[maxn], ft[maxn];
    	PPMdata **matrix;
    	int cnt = 0;
    	double w[maxn], Fw[maxn];
    	for (double w0 = -0.5 * PI - 0.5; w0 <= 0.5 * PI + 0.5; w0 += 0.01) {
    		w[cnt] = w0;
    		Fw[cnt] = F(w0);
    		++cnt;
    	}
    	arrayToPPM("out1.ppm", matrix, w, Fw, 1920, 1080, cnt, 0, 0, 3 * PI, 3, PI / 4.0, 0.5);
    	printf("------task 1 finished------\n");

    	cnt = 0;
    	for (double t0 = -20.0; t0 <= 20.0; t0 += 0.1) {
    		t[cnt] = t0;
    		ft[cnt] = fourierInvTransform(t0, F, -0.5 * PI - 0.5, 0.5 * PI + 0.5, 0.001);
    		++cnt;
    	}
    	arrayToPPM("out2.ppm", matrix, t, ft, 1920, 1080, cnt, 0, 0, 45, 1, 1, 0.03);
    	printf("------task 2 finished------\n");

    	double Fw_2[maxn], Fw_5[maxn], Fw_10[maxn];



    	cnt = 0;
    	for (double w0 = -10.0; w0 <= 10.0; w0 += 0.01) {
    		int j = 0;
    		for (double i = -20.0; i <= 20.0; i += 0.1) { //计算傅里叶积分
    			Fw_2[cnt] += cos(w0 * i) * ft[j] * sample(i, 0.2, 0.01) * 0.01;
    			Fw_5[cnt] += cos(w0 * i) * ft[j] * sample(i, 0.5, 0.01) * 0.01;
    			Fw_10[cnt] += cos(w0 * i) * ft[j] * sample(i, 1.0, 0.01) * 0.01;
    			++j;
    		}
    		w[cnt] = w0;
    		++cnt;
    	}
    	arrayToPPM("out3_1.ppm", matrix, w, Fw_2, 1920, 1080, cnt, 0, 0, 22, 0.1, PI / 4.0, 0.02);
    	arrayToPPM("out3_2.ppm", matrix, w, Fw_5, 1920, 1080, cnt, 0, 0, 22, 0.1, PI / 4.0, 0.02);
    	arrayToPPM("out3_3.ppm", matrix, w, Fw_10, 1920, 1080, cnt, 0, 0, 22, 0.1, PI / 4.0, 0.02);
    	printf("------task 3 finished------\n");

    	double Fw2_fil[maxn], Fw5_fil[maxn], Fw10_fil[maxn];
    	cnt = 0;
    	for (double w0 = -10.0; w0 <= 10.0; w0 += 0.01) {
    		Fw2_fil[cnt] = Fw_2[cnt] * filter(w0, 0.2, 0.01);
    		Fw5_fil[cnt] = Fw_5[cnt] * filter(w0, 0.5, 0.01);
    		Fw10_fil[cnt] = Fw_10[cnt] * filter(w0, 1.0, 0.01);
    		++cnt;
    	}
    	arrayToPPM("out4_1.ppm", matrix, w, Fw2_fil, 1920, 1080, cnt, 0, 0, 22, 4, PI / 4.0, 0.5);
    	arrayToPPM("out4_2.ppm", matrix, w, Fw5_fil, 1920, 1080, cnt, 0, 0, 22, 4, PI / 4.0, 0.5);
    	arrayToPPM("out4_3.ppm", matrix, w, Fw10_fil, 1920, 1080, cnt, 0, 0, 22, 4, PI / 4.0, 0.5);

    	cnt = 0;
    	double ft_2[maxn], ft_5[maxn], ft_10[maxn];
    	for (double t0 = -20.0; t0 <= 20.0; t0 += 0.1) {
    		int j = 0;
    		for (double w0 = -10.0; w0 <= 10.0; w0 += 0.01) {
    			ft_2[cnt] += cos(w0 * t0) * Fw2_fil[j] * 0.01;
    			ft_5[cnt] += cos(w0 * t0) * Fw5_fil[j] * 0.01;
    			ft_10[cnt] += cos(w0 * t0) * Fw10_fil[j] * 0.01;
    			++j;
    		}
    		ft_2[cnt] /= (2.0 * PI);
    		ft_5[cnt] /= (2.0 * PI);
    		ft_10[cnt] /= (2.0 * PI);
    		++cnt;
    	}
    	arrayToPPM("out4-4.ppm", matrix, t, ft_2, 1920, 1080, cnt, 0, 0, 45, 1, 1, 0.03);
    	arrayToPPM("out4-5.ppm", matrix, t, ft_5, 1920, 1080, cnt, 0, 0, 45, 1, 1, 0.03);
    	arrayToPPM("out4-6.ppm", matrix, t, ft_10, 1920, 1080, cnt, 0, 0, 45, 1, 1, 0.03);
    	return 0;
    }

  2. 其中的arratToPPM.h是本人编写的绘图用代码库。文件内容为：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    131
    132
    133
    134
    135
    136
    137
    138
    139
    140
    141
    142
    143
    144
    145
    146
    147
    148
    149
    150
    151
    152
    153
    154
    155
    156
    157
    158
    159
    160
    161
    162
    163
    164
    165
    166
    167
    168
    169
    170
    171
    172
    173
    174
    175
    176
    177
    178
    #include "stdio.h"
    #include "stdlib.h"
    #include "string.h"
    #include "math.h"

    #define INF 999999999
    #define linerFunc(x1,y1,x2,y2,x) ((x-x1)*(y2-y1)/(x2-x1)+y1)

    typedef struct {
    	int r, g, b;
    } PPMdata;

    int arrayToPPM(char *name,
                   PPMdata **matrix,
                   const double *x,
                   double *y,
                   int width,
                   int height,
                   int arrayLen,
                   double centerX,
                   double centerY,
                   double rangeX,
                   double rangeY,
                   double gridX,
                   double gridY);

    PPMdata makePPMdata(int r, int g, int b);
    int numToMatPos(double num, double center, double range, double picLen);
    void drawMatrix(PPMdata **matrix, int width, int height, int x, int y, PPMdata color);
    void drawPoint(PPMdata **matrix, int width, int height, int x, int y, PPMdata color, int size);
    void matToPPM(char *fileName, PPMdata **matrix, int width, int height);



    PPMdata makePPMdata(int r, int g, int b) {
    	PPMdata ans;
    	ans.r = r, ans.g = g, ans.b = b;
    	return ans;
    }

    int numToMatPos(double num, double center, double range, double picLen) {
    	return (int)(picLen / 2 + (num - center) / range * picLen);
    }

    void drawMatrix(PPMdata **matrix, int width, int height, int x, int y, PPMdata color) {
    	if (x >= width || x < 0 || y >= height || y < 0) return;
    	else matrix[y][x] = color;
    	return;
    }

    void drawPoint(PPMdata **matrix, int width, int height, int x, int y, PPMdata color, int size) {
    	for (int i = -1 * size; i <= size; ++i) {
    		for (int j = -1 * size; j <= size; ++j) {
    			int u = x + i, v = y + j;
    			if (u >= width || u < 0 || v >= height || v < 0) return;
    			else matrix[v][u] = color;
    		}
    	}
    }


    void matToPPM(char *fileName, PPMdata **matrix, int width, int height) {
    	freopen(fileName, "w", stdout);
    	printf("P3\n");
    	printf("%d %d\n", width, height);
    	printf("255\n");
    	for (int i = 0; i < height; ++i) {
    		for (int j = 0; j < width; ++j) {
    			printf("%d %d %d", matrix[i][j].r, matrix[i][j].g, matrix[i][j].b);
    			printf(" ");
    		}
    		printf("\n");
    	}
    	fclose(stdout);
    	freopen("CON", "w", stdout);
    	return;
    }

    int arrayToPPM(char *name,
                   PPMdata **matrix,
                   const double *x,
                   double *y,
                   int width,
                   int height,
                   int arrayLen,
                   double centerX,
                   double centerY,
                   double rangeX,
                   double rangeY,
                   double gridX,
                   double gridY) {

    	matrix = (PPMdata **) malloc(sizeof(PPMdata *) * height);
    	for (int i = 0; i < height; i++) {
    		matrix[i] = (PPMdata *) calloc(width, sizeof(PPMdata));
    	}
    	if (height % 2) ++height;
    	if (width % 2) ++width;
    	double *y_save=(double*)malloc(sizeof(double)*(arrayLen+2));
    	for(int i=0;i<arrayLen;++i) y_save[i]=y[i];
    	for (int i = 0; i < arrayLen; ++i) {
    		y[i] = centerY - (y[i] - centerY);
    	}
    	for (int i = 0; i < arrayLen; ++i) {
    		if (x[i] < x[i - 1] && i >= 1) {
    			fprintf(stderr, "ERROR:X is not increasing.\n");
    			return -1;
    		}
    		if (!(x[i] > centerX - (rangeX / 2.0) && x[i] < centerX + (rangeX / 2.0))) {
    			fprintf(stderr, "ERROR:X out of range.\n");
    			return -1;
    		}
    		if (!(y[i] > centerY - (rangeY / 2.0) && y[i] < centerY + (rangeY / 2.0))) {
    			if (y[i] > centerY + (rangeY / 2.0)) y[i] = centerY + (rangeY / 2.0);
    			if (y[i] < centerY - (rangeY / 2.0)) y[i] = centerY - (rangeY / 2.0);
    		}
    	}
    	PPMdata background, axis, grid, line;
    	background = makePPMdata(255, 255, 251);
    	axis = makePPMdata(28, 28, 28);
    	grid = makePPMdata(189, 192, 186);
    	line = makePPMdata(0, 92, 175);

    	for (int i = 0; i < height; ++i) {
    		for (int j = 0; j < width; ++j) {
    			matrix[i][j] = background;
    		}
    	}

    	for (int i = 0; centerX + i * gridX <= centerX + rangeX / 2; ++i) {
    		for (int j = 0; j < height; ++j) {
    			drawPoint(matrix, width, height, numToMatPos(centerX + i * gridX, centerX, rangeX, width), j, grid, 0);
    			drawPoint(matrix, width, height, numToMatPos(centerX - i * gridX, centerX, rangeX, width), j, grid, 0);
    		}
    	}
    	for (int i = 0; centerY + i * gridY <= centerY + rangeY / 2; ++i) {
    		for (int j = 0; j < width; ++j) {
    			drawPoint(matrix, width, height, j, numToMatPos(centerY + i * gridY, centerY, rangeY, height), grid, 0);
    			drawPoint(matrix, width, height, j, numToMatPos(centerY - i * gridY, centerY, rangeY, height), grid, 0);
    		}
    	}
    	for (int i = 0; i < width; ++i) {
    		drawPoint(matrix, width, height, i, height / 2, axis, 1);
    	}

    	for (int i = 0; i < height; ++i) {
    		drawPoint(matrix, width, height, width / 2, i, axis, 1);
    	}

    	double stepX = rangeX / width;
    	for (double X = centerX - rangeX / 2; X < centerX + rangeX / 2; X += stepX) {
    		int linerIndex = -1;
    		for (int i = 0; i < arrayLen; ++i) {
    			if (X < x[i]) {
    				linerIndex = i - 1;
    				break;
    			}
    		}
    		if (linerIndex == -1) {
    			continue;
    		} else {
    			double u, v, Y;
    			u = numToMatPos(X, centerX, rangeX, width);
    			Y = linerFunc(x[linerIndex], y[linerIndex], x[linerIndex + 1], y[linerIndex + 1], X);
    			v = numToMatPos(Y, centerY, rangeY, height);
    			drawPoint(matrix, width, height, u, v, line, 1);
    		}
    	}

    	matToPPM(name, matrix, width, height);
    	fprintf(stderr, "%s Draw finish!\n",name);
    	for(int i=0;i<height;++i)
    		free(matrix[i]);
    	free(matrix);
    	for(int i=0;i<arrayLen;++i) y[i]=y_save[i];
    	
    	return 0;
    }
实验
#信号与系统
信号与系统实验2：信号的矩形脉冲抽样与恢复
https://suzumiyaakizuki.github.io/2022/05/10/信号报告2/
作者
SuzumiyaAkizuki
发布于
2022年5月10日
许可协议