# 实验二：快速傅里叶变换

## 1. 给定序列，并对其进行相应的分析

### 1.1 高斯序列

x

a

(

n

)

=

{

e

(

n

p

)

2

q

0

n

15

0

x_a(n)=begin{cases} e^{-frac{(n-p)^2}{q}}~~ 0leq nleq 15 \ ~~ ~~ ~~0 ~~ ~~ ~~ ~~~~ ~~其他 \ end{cases}

xa(n)={eq(np)2  0n15      0

clc;clear;
n = 0:15; p1 = 8; p2 = 13; p3 = 14; q1 = 2; q2 = 4; q3 = 8;
% p = 8, q = 2, 4, 8
x1 = exp(-(n-p1).*(n-p1)/q1); x2 = exp(-(n-p1).*(n-p1)/q2); x3 = exp(-(n-p1).*(n-p1)/q3);
x1w = fft(x1); x2w = fft(x2); x3w = fft(x3);
figure(1);
subplot(3,2,1); stem(n, x1, 'fill'); ylabel('p=8, q=2'); title('时域');
subplot(3,2,2); stem(n, abs(x1w), 'fill'); title('频域');
subplot(3,2,3); stem(n, x2, 'fill'); ylabel('p=8, q=4')
subplot(3,2,4); stem(n, abs(x2w), 'fill');
subplot(3,2,5); stem(n, x3, 'fill'); ylabel('p=8, q=6')
subplot(3,2,6); stem(n, abs(x3w), 'fill');
sgtitle('高斯序列时域和幅频特性');
% q = 8, p = 8, 13, 14
x1 = exp(-(n-p1).*(n-p1)/q3); x2 = exp(-(n-p2).*(n-p2)/q3); x3 = exp(-(n-p3).*(n-p3)/q3);
x1w = fft(x1); x2w = fft(x2); x3w = fft(x3);
figure(2);
subplot(3,2,1); stem(n, x1, 'fill'); ylabel('p=8, q=8'); title('时域');
subplot(3,2,2); stem(n, abs(x1w), 'fill'); title('频域');
subplot(3,2,3); stem(n, x2, 'fill'); ylabel('p=13, q=8')
subplot(3,2,4); stem(n, abs(x2w), 'fill');
subplot(3,2,5); stem(n, x3, 'fill'); ylabel('p=14, q=8')
subplot(3,2,6); stem(n, abs(x3w), 'fill');
sgtitle('高斯序列时域和幅频特性');


### 1.2 衰减正弦序列

x

b

(

n

)

=

{

e

a

n

s

i

n

(

2

π

f

n

)

0

n

15

0

x_b(n)=begin{cases} e^{-an}sin(2pi fn)~~ 0leq nleq 15 \ ~~ ~~ ~~ ~~0 ~~ ~~ ~~ ~~ ~~其他 \ end{cases}

xb(n)={eansin(2πfn)  0n15        0

n=0:1:15; a=0.1; f1=0.0625; f2=0.04375; f3=0.05625;
xb1=exp(-a*n).*sin(2*pi*f1*n);
figure
subplot(3, 2, 1); stem(n, xb1,'.');
title("f=0. 0625的时域特性"); xlabel('n'); ylabel('xb1 (n)'); grid on;
[H, w] = freqz(xb1,  1,  [], 'whole', 1);
Hamplitude = abs (H);
subplot(3, 2, 2); plot (w, Hamplitude);
title("f=0. 0625的幅频响应"); xlabel('w/ (2*pi)'); ylabel('|H(exp(jw))|');
grid on
xb2=exp(-a*n) .* sin(2*pi*f2*n) ;
subplot(3, 2, 3); stem(n, xb2, '.');
title("f=0. 04375的时域特性"); xlabel('n'); ylabel('xb2 (n)'); grid on
[H, w] = freqz(xb2, 1, [], 'whole', 1);
Hamplitude = abs (H);
subplot(3, 2, 4);
plot(w, Hamplitude);
title("f=0. 04375的幅频响应"); xlabel('w/(2*pi)' ); ylabel('|H(exp(jw))|'); grid on
xb3=exp(-a*n) .* sin (2*pi*f3*n);
subplot (3, 2, 5);
stem(n, xb3,'.' );
title("f=0. 05625的时域特性"); xlabel('n'); ylabel('xb3(n)');
grid on
[H, w] = freqz(xb3, 1, [], 'whole', 1);
Hamplitude = abs (H);
subplot(3, 2, 6); plot (w,Hamplitude);
title("f=0. 05625的幅频响应"); xlabel('w/(2*pi)'); ylabel('|H(exp(jw))|'); grid on


f=0.0625时，谱峰位置出现正确，但是采样频率 fs = 1/T 小于两倍的信号频率，存在在混叠现象，时域采样为一周期，不满足采样定理。

### 1.3 三角波序列

x

c

(

n

)

=

{

n

0

n

3

8

n

4

n

7

0

x_c(n)=begin{cases} n~~~~~~~~~~~0leq nleq 3 \ 8-n~~~~4leq nleq 7 \ 0~~~~~~~~~~~其他 \ end{cases}

xc(n)=n           0n38n    4n70

### 1.4 2反三角波序列

x

d

(

n

)

=

{

4

n

0

n

3

n

4

4

n

7

0

x_d(n)=begin{cases} 4-n~~~~0leq n leq 3 \ n-4~~~~4leq n leq 7 \ 0~~~~~~~~~~~其他 end{cases}

xd(n)=4n    0n3n4    4n70

close all; clear; clc;
n1 = 0:1:3; xc1 = n1+1; n2 = 4:7; xc2 = 8-n2;
xc = [xc1,xc2]; n = [n1,n2];
subplot(321)
stem(n,xc); xlabel ('n'); ylabel('xc'); title('三角序列');

n1 = 0:1:3; xd1 = 4-n1; n2 = 4:7; xd2 = n2-3;
xd = [xd1, xd2]; n = [n1,n2];
subplot(322)
stem(n,xd); xlabel('n') ; ylabel ('xd') ; title('反三角序列');

N=16;
[H1,w1] = freqz (xc,1, 256, 'whole', 1);
Hamplitude1 = abs (H1);
subplot(323)
plot (2*w1, Hamplitude1);
title('xc幅频响应'); xlabel ('w/pi'); ylabel ('IH(exp(jw)) 1'); grid on

[H2,w2] = freqz (xd,1, 256, 'whole', 1) ;
Hamplitude2 = abs (H2);
subplot(324)
plot(2*w2, Hamplitude2);
title('xd幅频响应'); xlabel ('w/pi'); ylabel('|H(exp(jw))|'); grid on

[H3, w3] = freqz(xc, 1, N, 'whole', 1);
Hamplitude3 = abs (H3) ;
subplot(325);
h3 = stem(2*w3, Hamplitude3, '*');
title ( 'xc幅频响应进行N点FFT'); xlabel ('n'); ylabel(' |H(exp(jw))|'); grid on

[H4, w4] = freqz(xd, 1, N, 'whole', 1) ;
Hamplitude4 = abs (H4);
subplot(326);
h4 = stem(2*w4, Hamplitude4, '*');
title ( 'xd幅频响应进行N点FFT'); xlabel ('n'); ylabel('|H(exp(jw))|'); grid on


## 2. 分析有限长正弦序列

clear all; clc; close all;
f1=10; f2=11; fs=64; N=32;
dt=1/fs;
T=(0:N-1)*dt;
x1=sin(2*pi*f1*T);
x2=sin(2*pi*f2*T);
x1_dft=fft(x1); x2_dft=fft(x2);
subplot(221); plot(T,abs(x1_dft/N), '-*'); xlabel('频率/HZ'); ylabel('振幅'); title('f=10HZ DFT变换');
subplot(222); plot(T,abs(x2_dft/N), '-*'); xlabel('频率/HZ'); ylabel('振幅'); title('f=11HZ DFT变换');

[h1,w1]=freqz(x1,1); [h2,w2]=freqz(x2,7);
subplot(223); plot(w1,abs(h1)); xlabel('频率/HZ'); ylabel('振幅'); title('f=10HZ DTFT变换');
subplot(224); plot(w2,abs(h2)); xlabel('频率/HZ'); ylabel('振幅'); title('f=11HZ DTFT变换');


## 实验要求：

1. 熟悉MATLAB的信号处理函数
2. 熟练离散信号时域和频域图形的绘制
3. 熟练绘制信号的幅频响应和相频响应
4. 完成频谱与实例分析

## 附：相关MATLAB函数

y =exp( X)；对向量X的各元素作指数计算，结果为一向量。

conj(X)；对向量X的各元素作复共辄运算，即将虚部改变符号。

real(X)；对向量X的各元素取其实部。

v =randn ( size(X))；生成同X具有相同维数的正态(高斯)分布的随机矩阵，通常用于生成高斯白噪声。

fft(X，N)；计算X的N点FFT，如果X的长度小于N，则将在X后补零。反之，如果X的长度大于N，则将对X进行截取。

ifft(X)；计算X的N点 IFFT。

fftshift(Y)；如果Y为向量，该命令将Y分成左右两部分并交换位置。

xcorr(A ，B)；计算A、B两个序列的互相关。
n ( size(X))；生成同X具有相同维数的正态(高斯)分布的随机矩阵，通常用于生成高斯白噪声。

fft(X，N)；计算X的N点FFT，如果X的长度小于N，则将在X后补零。反之，如果X的长度大于N，则将对X进行截取。

ifft(X)；计算X的N点 IFFT。

fftshift(Y)；如果Y为向量，该命令将Y分成左右两部分并交换位置。

xcorr(A ，B)；计算A、B两个序列的互相关。

THE END