吉布斯现象作业

Posted on | In

使用matlab和C进行绘图与计算。

程序如下：
C程序计算采样和DFT频域转换

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
#include <stdio.h>
#include <math.h>
//# define NP 2001
# define PI 2*asin(1.0)
#define T 2.0
void dft(xreal, ximag);
float xtime(float dt, float *x, float *t, int N);

int main()
{
    float dt1=1.6;
    float dt2=2.5;
    int N1=T*1000/dt1;
    printf("%d ",N1);

    float x1[N1],t1[N1];
    int N2 = T*1000/dt2;
    float x2[N2],t2[N2];

    xtime(dt1,x1,t1, N1);
    xtime(dt2,x2,t2, N2);

    float xreal1[N1], ximag1[N1];
    printf("输出测试");  //测试无误
    for(int i=0; i<N1; i++)
    {
        xreal1[i]=x1[i]; ximag1[i]=0;  //将雷克子波输入实部为原值，虚部为0
        printf("xreal1 %f\t%d\n", xreal1[i], i);
    }
    dft(xreal1, ximag1,N1);
    for(int i=0; i<N1; i++)
    {
        printf("%f\t%f\t%d\n", xreal1[i],ximag1[i],i);
    }

    float xreal2[N1], ximag2[N1];
    printf("输出测试");  //测试无误
    for(int i=0; i<N2; i++)
    {
        xreal2[i]=x1[i]; ximag2[i]=0;  //将雷克子波输入实部为原值，虚部为0
        printf("xreal2 %f\t%d\n", xreal2[i], i);
    }
    dft(xreal2, ximag2, N2);
    for(int i=0; i<N1; i++)
    {
        printf("%f\t%f\t%d\n", xreal2[i],ximag2[i], i);
    }

    FILE *DFT;  //打开最终输出的文件
    DFT = fopen("work_dt1.csv", "w");
printf( "dt1 \n");
    for(int i=0; i<N1; i++)
    {
        printf("%10.6f\t %10.6f\n", sqrt(pow(xreal1[i],2)+pow(ximag1[i], 2)), i/dt1/N1);
         fprintf(DFT, "%10.6f\t %10.6f\n", sqrt(pow(xreal1[i],2)+pow(ximag1[i], 2)), i/dt1/N1);
    }
    FILE *DFT2;
    DFT2 = fopen("work_dt2.csv", "w");
    printf( "dt2 \n");
     for(int i=0; i<N2; i++)
    {
         printf("%10.6f\t %10.6f\n", sqrt(pow(xreal2[i],2)+pow(ximag2[i], 2)), i/dt2/N2);
         fprintf(DFT2, "%10.6f\t %10.6f\n", sqrt(pow(xreal2[i],2)+pow(ximag2[i], 2)), i/dt2/N2);
    }
fclose(DFT);
fclose(DFT2);

    return 0;
}

void dft(float *xreal,float *ximag, int N)  //DFT 离散傅里叶正变换
{
    float Xreal[N],Ximag[N];
    float sumreal, sumimag, rad;  //定义实部计算加和、虚部计算加和、角度

    for(int k = 0; k < N; k++)
    {
        sumreal = 0.0; sumimag = 0.0;
        for(int n=0; n<N; n++)
        {
            rad = 2*PI/N*n*k;
            //printf("thisxreal %f\t%d\t%d\n", xreal[n], n, k);
            sumreal += xreal[n]*cos(rad) + ximag[n]*sin(rad);
            sumimag += -xreal[n]*sin(rad) + ximag[n]*cos(rad);

        }

        Xreal[k] = sumreal;
        Ximag[k] = sumimag;
    }
    for(int k = 0; k < N; k++)
    {
        xreal[k] = Xreal[k];
        ximag[k] = Ximag[k];
    }
}

float xtime(float dt, float *x, float *t, int N)
{

    for(int i=0;i<N;i++)
    {
        t[i] = -T/2+i*dt/1000;
        if(t[i]==0)
             x[i] = 227.0;
        else
        x[i] = sin(300*t[i])/PI/t[i];
    }

}

MATLAB进行图像绘制

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
dt1 = 1.6;
dt2 = 2.5;
[t1,x1,N1] = xtime(dt1);
[t2,x2,N2] = xtime(dt2);

figure(1)
hold on
subplot(2,1,1)
plot(t1,x1,'-r')
legend('采样间隔1.6ms');
title('时域采样结果')
subplot(2,1,2)
plot(t2,x2,'-b')
legend('采样间隔2.5ms');

load work_dt1.csv
load work_dt2.csv

figure(2)
subplot(2,1,1)
hold on
plot(work_dt1(:,2),work_dt1(:,1),'-r')
legend('采样间隔1.6ms')
title('频谱绘制')
subplot(2,1,2)
hold on
plot(work_dt2(:,2),work_dt2(:,1),'-b')
legend('采样间隔2.5ms')

function [t,x,N] = xtime(dt)
    T = 2; % 采样总时间为2s，原点两边各采样1s
    N = floor(T*1000/dt);
    x = zeros(1,N);
    t = zeros(1,N);
    for i=1 : N
        t(i) = -T/2+(i-1)*dt/1000;
        x(i) = 1/pi/t(i)*sin(300*t(i));
    end
end

频率采样

时域采样