实验一：数字信号的产生和DTFT运算

因为现实世界里存在的是模拟信号，因此数字信号处理的第一个问题是将信号离散化，得到一个数字信号，然后再进行数字处理。

(1) 常用数字信号序列的产生：

熟悉 Matlab 产生数字信号的基本命令，加深对数字信号概念的理解，并能够用 Matlab 产生和绘制出一些常用离散信号序列，例如 δ(n)、单位阶跃序列 u(n)、矩形序列 R(n)、正弦序列 Sin(nw) 等。

(2) 数字信号的基本运算：

加、减、尺度（乘除）和移位是数字信号处理中最基本的算术运算，将上述基本序列进行这些基本运算，得到多个序列构成的组合序列。

通过本次实验，掌握 Matlab 中这些基本运算命令，对数字信号处理有一个基本概念，为后面的数字信号分析和滤波打下基础。

数字信号处理的一个重要分支就是信号分析，而信号分析的基本工具是离散傅立叶变换。利用傅立叶变换和级数所形成的频谱分析技术作为处理连续信号的重要工具已经应用得很久了,1956年库力（Cooley）和图基（Tukey）所发展的近似频谱的快速算法为频谱分析的数字信号的谱分析铺平了道路。因此，DFT（FFT）得到广泛应用。

（3）假设信号 x(n) 由下述信号组成：

$x(n)=0.1*cos(0.45n\pi)+sin(0.3n\pi)-cos(0.31n\pi=\frac{\pi}{4})$

这个信号有两根主谱线 0.3pi 和 0.31pi 靠的非常近，而另一根谱线 0.45pi 的幅度很小，请选择合适的长度 N 和窗函数，用 DTFT 和 DFT 分析其频谱，得到清楚的三根谱线。

通过本次实验，应该掌握：

(a) 用傅立叶变换进行信号分析时基本参数的选择。

(b) 经过离散时间傅立叶变换（DTFT）和有限长度离散傅立叶变换（DFT）后信号频谱上的区别，前者 DTFT 时间域是离散信号，频率域还是连续的，而 DFT 在两个域中都是离散的。

%% FILE INTRODUCTION 
% file name:    sequenceGenerate.m
% funtion:      Digital Signal Processing Experiment 1(1)
% author:       chengkai
% date:         2016/12/1
% release:      First Release

clc,clear % clear the command window and clear variables in the workspace window
close all % close all figure widows
% 产生单位采样序列
delta_sequence = -10:10;
delta_n = delta_sequence == 0;
stem(delta_sequence, delta_n);
axis([-12 12 -0.1 1.5]) % 设置坐标轴
xlabel('n')
ylabel('\delta(n)')
title('单位采样序列\delta(n)')

% 产生单位阶跃序列u(n)
un_sequence = -10:10;
un = (un_sequence >= 0);
figure
stem(un_sequence, un)
axis([-12 12 -0.1 1.5]) % 设置坐标轴
xlabel('n')
ylabel('u(n)')
title('单位阶跃序列u(n)')

% 产生矩形序列R(n)
Rn_sequence = -10:10;
Rn = (Rn_sequence >= 0 & Rn_sequence <= 3); % 注意：此处必须是&而不能是&&
figure
stem(Rn_sequence, Rn)
axis([-12 12 -0.1 1.5]) % 设置坐标轴
xlabel('n')
ylabel('R_4(n)')
title('矩形序列R_{4}(n)')

% 产生正弦序列sin(nw)
sin_sequence = 0:99;
w = 0.04 * pi;
sin_nw = sin(sin_sequence * w);
figure
stem(sin_sequence, sin_nw)
axis([0 99 -1.5 1.5]) % 设置坐标轴
xlabel('n')
ylabel('sin(0.5\pi n)')
title('正弦序列sin(0.5\pi n)')

%% FILE INTRODUCTION

% file name: sequenceGenerate.m

% funtion: Digital Signal Processing Experiment 1(1)

% author: chengkai

% date: 2016/12/1

% release: First Release

clc,clear % clear the command window and clear variables in the workspace window

close all % close all figure widows

% 产生单位采样序列

delta_sequence = -10:10;

delta_n = delta_sequence == 0;

stem(delta_sequence, delta_n);

axis([-12 12 -0.1 1.5]) % 设置坐标轴

xlabel('n')

ylabel('\delta(n)')

title('单位采样序列\delta(n)')

% 产生单位阶跃序列u(n)

un_sequence = -10:10;

un = (un_sequence >= 0);

figure

stem(un_sequence, un)

axis([-12 12 -0.1 1.5]) % 设置坐标轴

xlabel('n')

ylabel('u(n)')

title('单位阶跃序列u(n)')

% 产生矩形序列R(n)

Rn_sequence = -10:10;

Rn = (Rn_sequence >= 0 & Rn_sequence <= 3); % 注意：此处必须是&而不能是&&

figure

stem(Rn_sequence, Rn)

axis([-12 12 -0.1 1.5]) % 设置坐标轴

xlabel('n')

ylabel('R_4(n)')

title('矩形序列R_{4}(n)')

% 产生正弦序列sin(nw)

sin_sequence = 0:99;

w = 0.04 * pi;

sin_nw = sin(sin_sequence * w);

figure

stem(sin_sequence, sin_nw)

axis([0 99 -1.5 1.5]) % 设置坐标轴

xlabel('n')

ylabel('sin(0.5\pi n)')

title('正弦序列sin(0.5\pi n)')

sequencegenerate_01 sequencegenerate_02

sequencegenerate_03

sequencegenerate_04

%% FILE INTRODUCTION 
% file name:    lab1_DTFT_DFT.m
% funtion:      Digital Signal Processing Experiment 1(3)
% author:       chengkai
% date:         2016/12/1
% release:      First Release

clc,clear % clear the command window and clear variables in the workspace window
close all % close all figure widows
N = 400; % 200为最小周期
n = 0:N-1;
w = 0 : 0.001 : 2 * pi; %待观察的频率向量
xn = 0.1 * cos(0.45 * n * pi) + sin(0.3 * n * pi) - cos(0.31 * n * pi - pi / 4);
subplot(311)
stem(n, xn)

tic
Xejw = DTFT(xn, w);
toc
subplot(312)
plot(w, abs(Xejw))
axis([0 max(w) 0 max(abs(Xejw) + 50)])

[Xk, k] = DFT(xn, N);
subplot(313)
stem(k,abs(Xk))
axis([0 max(k) 0 max(abs(Xk) + 50)])

% tic
% Xejw2 = DTFT_slow(xn, w);
% toc
% figure
% plot(w, abs(Xejw2))
% error = abs(Xejw2) - abs(Xejw);

%% 高密度谱与高分辨率谱的区别解析
% 补零并不能提高DFT的分辨率，补零可以得到在当前分辨率下的高密度频谱
% DFT的频谱分辨率是指对信号中两个靠的比较近的频谱分量的识别能力，它仅取决于截取连续信号的长度，在采样
% 频率不变的前提下，可以通过改变采样点数N来改变DFT的分辨率。
% 高密度频谱是指当信号的时间与长度不变时，在频域内对它的频谱提高采样率而得到的高密度谱，克服栅栏效应，但
% 不能提高DFT的分辨率。
figure
xnn = [xn zeros(1,N)];
Xejw2 = DTFT(xnn, w);
subplot(211)
plot(w, abs(Xejw2))
axis([0 max(w) 0 max(abs(Xejw2) + 50)])

subplot(212)
Xk2 = DFT(xnn, 2 * N);
subplot(313)
stem(0:2 * N-1,abs(Xk2))

%% FILE INTRODUCTION

% file name: lab1_DTFT_DFT.m

% funtion: Digital Signal Processing Experiment 1(3)

% author: chengkai

% date: 2016/12/1

% release: First Release

clc,clear % clear the command window and clear variables in the workspace window

close all % close all figure widows

N = 400; % 200为最小周期

n = 0:N-1;

w = 0 : 0.001 : 2 * pi; %待观察的频率向量

xn = 0.1 * cos(0.45 * n * pi) + sin(0.3 * n * pi) - cos(0.31 * n * pi - pi / 4);

subplot(311)

stem(n, xn)

tic

Xejw = DTFT(xn, w);

toc

subplot(312)

plot(w, abs(Xejw))

axis([0 max(w) 0 max(abs(Xejw) + 50)])

[Xk, k] = DFT(xn, N);

subplot(313)

stem(k,abs(Xk))

axis([0 max(k) 0 max(abs(Xk) + 50)])

% tic

% Xejw2 = DTFT_slow(xn, w);

% toc

% figure

% plot(w, abs(Xejw2))

% error = abs(Xejw2) - abs(Xejw);

%% 高密度谱与高分辨率谱的区别解析

% 补零并不能提高DFT的分辨率，补零可以得到在当前分辨率下的高密度频谱

% DFT的频谱分辨率是指对信号中两个靠的比较近的频谱分量的识别能力，它仅取决于截取连续信号的长度，在采样

% 频率不变的前提下，可以通过改变采样点数N来改变DFT的分辨率。

% 高密度频谱是指当信号的时间与长度不变时，在频域内对它的频谱提高采样率而得到的高密度谱，克服栅栏效应，但

% 不能提高DFT的分辨率。

figure

xnn = [xn zeros(1,N)];

Xejw2 = DTFT(xnn, w);

subplot(211)

plot(w, abs(Xejw2))

axis([0 max(w) 0 max(abs(Xejw2) + 50)])

subplot(212)

Xk2 = DFT(xnn, 2 * N);

subplot(313)

stem(0:2 * N-1,abs(Xk2))

%% FILE INTRODUCTION 
% file name:    DTFT.m
% funtion:      Caculate the DTFT of a discrete sequence using matrix
%               caculation to reduce time complexity.
% author:       chengkai
% date:         2016/12/1
% release:      First Release

function [ Xejw ] = DTFT( xn,w )
%DTFT 计算离散时间傅里叶变换
%   xn：待分析的序列
%   w:  待观察的频率向量
    xn_cnt=length(xn);%序列的长度
    n=0:xn_cnt-1;
    Xejw=xn*(exp(-1i).^(n'*w));
end

%% FILE INTRODUCTION

% file name: DTFT.m

% funtion: Caculate the DTFT of a discrete sequence using matrix

% caculation to reduce time complexity.

% author: chengkai

% date: 2016/12/1

% release: First Release

function [ Xejw ] = DTFT( xn,w )

%DTFT 计算离散时间傅里叶变换

% xn：待分析的序列

% w: 待观察的频率向量

xn_cnt=length(xn);%序列的长度

n=0:xn_cnt-1;

Xejw=xn*(exp(-1i).^(n'*w));

end

%% FILE INTRODUCTION 
% file name:    DFT.m
% funtion:      Caculate the DFT of a discrete sequence using matrix
%               caculation to reduce time complexity.
% author:       chengkai
% date:         2016/12/1
% release:      First Release
% other:        If the length of xn is less than N, the program will
%               add zeros behind xn sequence.

function [ Xk,k ] = DFT( xn,N )
%DFT 计算序列的N点DFT
%   xn：待分析的序列
%   N：采样点数
%   k:为归一化频率
    k = 0 : N - 1;
    n = 0 : N - 1;
    xn_cnt = length(xn);
    if xn_cnt < N
        yn = zeros(1,N);
        yn(1, 1 : xn_cnt) = xn;
        xn = yn;
    end
    Xk = xn(1, 1 : N) * (exp(-1i) .^ (2 * pi * n' * k / N));
    k = 2 * pi / N * k;
end

%% FILE INTRODUCTION

% file name: DFT.m

% funtion: Caculate the DFT of a discrete sequence using matrix

% caculation to reduce time complexity.

% author: chengkai

% date: 2016/12/1

% release: First Release

% other: If the length of xn is less than N, the program will

% add zeros behind xn sequence.

function [ Xk,k ] = DFT( xn,N )

%DFT 计算序列的N点DFT

% xn：待分析的序列

% N：采样点数

% k:为归一化频率

k = 0 : N - 1;

n = 0 : N - 1;

xn_cnt = length(xn);

if xn_cnt < N

yn = zeros(1,N);

yn(1, 1 : xn_cnt) = xn;

xn = yn;

end

Xk = xn(1, 1 : N) * (exp(-1i) .^ (2 * pi * n' * k / N));

k = 2 * pi / N * k;

end

lab1_dtft_dft_01

lab1_dtft_dft_02

分享到：