欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识！
您现在的位置是： 首页

【FPGA作业】第二章

程序员文章站 2022-07-14 20:02:19
...

二 正弦信号频谱分析实验

2.1 单音正弦信号采样序列的时域绘图

绘制下列信号的时域图

信号采样率fs=8000Hz
信号采样序列长度N=32
配置参数
f1=1000; Amp1=1; phy1=0
f2=7000; Amp2=1; phy2=0
f3=9000; Amp3=1; phy3=0

2.1.1 matlab程序

主程序，包括参数设置lab1.m

% lab1
close all;
clear;
fig_fname = 'sine_1_tone.jpg';
fs = 8E3; % 采样率 8k
N = 32; %向量长度
% frequency & amp & initial phase
f1 = 1000; Amp1 = 1; phy1 = 0;
f2 = 7000; Amp2 = 1; phy2 = 0;
f3 = 9000; Amp3 = 1; phy3 = 0;
% 信号向量的下标索引
n_idx = [0 : N-1];
% generate signal 
x1 = Amp1 * sin( 2*pi*f1/fs*n_idx + phy1);
x2 = Amp2 * sin( 2*pi*f2/fs*n_idx + phy2);
x3 = Amp3 * sin( 2*pi*f3/fs*n_idx + phy3);
% 生成绘图的标题文字 
str_fs = num2str(fs); str_N = num2str(N);
str_f1 = num2str(f1); str_Amp1 = num2str(Amp1); str_phy1 = num2str(phy1);
str_f2 = num2str(f2); str_Amp2 = num2str(Amp2); str_phy2 = num2str(phy2);
str_f3 = num2str(f3); str_Amp3 = num2str(Amp3); str_phy3 = num2str(phy3);
% 合成绘图的标题文字
title_str1 = ['f1:', str_f1, 'Hz, ','fs:', str_fs, 'Hz, ','N:', str_N, ', Amp1:' str_Amp1, ', phy1:',str_phy1];
title_str2 = ['f2:', str_f2, 'Hz, ','fs:', str_fs, 'Hz, ','N:', str_N, ', Amp2:' str_Amp2, ', phy2:',str_phy2];
title_str3 = ['f3:', str_f3, 'Hz, ','fs:', str_fs, 'Hz, ','N:', str_N, ', Amp2:' str_Amp3, ', phy2:',str_phy3];
%绘图
h = figure;

subplot(3,1,1);
plot(n_idx, x1, 'color', 'blue'); %连线方式绘图
grid on; hold;
stem(n_idx, x1, 'color', 'red'); %离散样值方式绘图
title(title_str1, 'fontsize', 14);

subplot(3,1,2);
plot(n_idx, x2, 'color', 'blue'); %连线方式绘图
grid on; hold;
stem(n_idx, x2, 'color', 'red'); %离散样值方式绘图
title(title_str2, 'fontsize', 14);

subplot(3,1,3);
plot(n_idx, x3, 'color', 'blue'); %连线方式绘图
grid on; hold;
stem(n_idx, x3, 'color', 'red'); %离散样值方式绘图
title(title_str3, 'fontsize', 14);

2.1.2 三种信号时域图

【FPGA作业】第二章

图2. 1 1K、7K、9KHz信号时域图

2.1.3 实验结果分析

根据采样定理，信号的频率为f0，采样频率fs，当fs>2f0时，可将信号完全重构。本实验中三种频率的单音信号1kHz、7kHz和9kHz，采样频率为8kHz，由图看出，时序图相似。分析为三种信号频率差异过小，当频率为500Hz的信号被8kHz的采样率重构时，可看出采样率对不同频率信号的影响。
【FPGA作业】第二章
图2. 2 500、7K、9KHz信号时域图

2.2 单音正弦信号采样序列的DFT结果绘图

2.2.1 实验要求

-计算已知正弦信号的DFT谱
-设定不同的仿真配置，观察正弦信号的DFT谱
* 改变DFT的计算长度
* 改变输入正弦信号采样序列的频率

分析信号如下
采样率fs=8kHz
信号采样序列长度N=32
配置参数
f1=3000 Hz; Amp1=1; phy1=0
f2=500Hz; Amp2=1; phy2=0
f3=3000Hz; Amp3=1; phy3=0

2.2.2 matlab程序

主程序，包括参数设置sine_1_tone_dft.m

% file: sine_1_tone_dft.m
% test with Matlab r2015b
close all; clear;
fig_fname_0 = 'sine_1_tone_dft_0.jpg';
fig_fname_1 = 'sine_1_tone_dft_1.jpg';
fig_fname_2 = 'sine_1_tone_dft_2.jpg';
fs  = 8E3; 
f0_0  = 3000; f0_1  = 500; f0_2  = 3000; 
N_0   = 32 ; N_1   = 32 ; N_2   = 35 ;

% plot each case
func_sine_1_tone_dft_plot(f0_0, fs, N_0, fig_fname_0);
func_sine_1_tone_dft_plot(f0_1, fs, N_1, fig_fname_1);
func_sine_1_tone_dft_plot(f0_2, fs, N_2, fig_fname_2);
绘制信号的时域图、线性幅度谱和对数幅度谱func_sine_1_tone_dft_plot.m
function func_sine_1_tone_dft_plot(f0, fs, N, fig_fname);
  n_idx = [0:N-1];  % 信号序列的下标向量 
  x = sin(2*pi*f0/fs*n_idx);    % 生成正弦信号序列
  y = abs(fft(x)); y_dB = 20*log10(1E-6+y); % 计算 DFT 结果线性幅度谱及对数幅度谱

  h = figure;
  % 绘图横轴左右边界
  x_left = -1; x_right = N+1;
  % 时域绘图 
  subplot(3,1,1);
  stem(n_idx, x); grid on; hold;
  plot(n_idx, x); xlim([x_left, x_right]);
  xlim([x_left, x_right]);
  title_str = ['f0:', num2str(f0),'Hz, ', 'fs:', num2str(fs),'Hz, ' 'N:', num2str(N)];
  title(title_str, 'fontsize',14);
  % 绘图线性幅度谱
  subplot(3,1,2); stem(n_idx, y); grid on;
  xlim([x_left, x_right]);
  title_str = ['DFT Magtitude Linear Scale, '];
  title_str = [title_str, ' Max:', num2str(max(y)), ', Min:', num2str(min(y))];
  title(title_str, 'fontsize',14);
  % 绘图对数幅度谱
  subplot(3,1,3); stem(n_idx, y_dB);  grid on;
  xlim([x_left, x_right]);
  title_str = ['DFT Magtitude dB Scale, '];
  title_str = [title_str, ' Max:', num2str(max(y_dB)), ', Min:', num2str(min(y_dB))];
  title(title_str, 'fontsize',14);

  print(h, '-djpeg', fig_fname);  % 保存绘图结果文件
end

2.2.3 实验结果

【FPGA作业】第二章
图2. 3 3KHz时域、频域、对数频谱图
采样序列中包含12个完整的信号周期采样值，DFT幅度谱上有两根峰值谱线，其余位置均为0.
【FPGA作业】第二章
图2. 4 500Hz信号时域、频域、对数频谱图

采样序列包括2个完整的采样周期，DFT幅度谱有两根峰值谱线，其余位置均为0.
【FPGA作业】第二章

图2. 5 3KHz信号时域、频域、对数频谱图
采样序列中包含13个完整的信号周期，以及一个不完整的信号周期采样值。DFT幅度谱分析结果有2根极大值谱线，但其余位置不为0.

2.2.4 实验分析

当采样周期是信号周期的整数倍时，信号采样后的DFT分析结果与连续信号的傅里叶变换频谱最相似。否则发生频谱泄露。

2.3 有限长矩形窗序列的DTFT

因计算机只能计算有限长和离散化的数据，因此使用计算机进行信号处理时，需要将模拟信号变为有限长、离散的数字信号，即对信号进行离散事件序列傅里叶变换DTFT。

2.3.1 matlab程序

主程序，包括参数设置rec_dtft_plot.m

clc 
clear
close all
w_min = -1*pi;
w_max =  1*pi;
w_delta = pi/1000;
% 可修改
n_min =   0;
n_max =   3 ;

n = [n_min:1:n_max];
w = [w_min:w_delta:w_max];

N = length(n);

x = ones(1,N); 
x_dscp_str = 'x(n)';

Y = func_calc_dtft(x, n, w); 

h = func_plot_dtft(x,n,w,Y, x_dscp_str);
计算矩形窗的DTFT, func_calc_dtft.m
% Y(w) = DTFT(x[n]);
function Y = func_calc_dtft(x, n, w) 

N     = length(n);
N_w   = length(w);

Y = zeros(N_w,1);
for idx = 1:N_w
  e_jw = exp(-j*w(idx)*n);
  Y(idx) = sum(e_jw .* x);
end

绘图，包括矩形窗函数的时域图、频谱实部、频谱虚部和频谱的模值

func_plot_dtft.m
% Y[w] = DTFT(x[n]);
% plot time sequence x[n], real/imag part of Y[w], abs of Y[w]
function h = func_plot_dtft(x,n,w,Y, x_dscp_str);
h = figure; 

N = length(n);
n_min = min(n);
n_max = max(n);
w_min = min(w);
w_max = max(w);
N_w = length(w);
w_delta = (w_max-w_min)/(N_w-1);

w_x_tick_delta = 1*pi/4;

w_x_tick = get_w_x_tick_vec(w_min, w_max, w_x_tick_delta);

w_x_tick_label_cell = get_w_x_tick_label_cell(w_x_tick, w_x_tick_delta);

Y_real = real(Y); Y_imag = imag(Y); Y_abs = abs(Y);

subplot(2,2,1);
title_str = [x_dscp_str, ', n \in [',int2str(n_min),',', int2str(n_max), ... 
            '], ', int2str(N), ' samples'];
stem(n, x, 'fill'); title(title_str, 'fontsize',14);
y_lim_max = 1.1*max(abs(x)); y_lim_min = -y_lim_max; 
grid on; xlim([n_min-1,n_max+1]); ylim([y_lim_min, y_lim_max]);

subplot(2,2,2);
Left_str  = get_tick_str((w_min/pi),1,'\pi');
Right_str = get_tick_str((w_max/pi),1,'\pi');
w_xLabelStr =['\omega: [', Left_str,',', Right_str, ']'];
title_str = ['\{X(e^{j\omega})\}, \Delta\omega = \pi/', num2str(pi/w_delta)];
plot(w, Y_real);title(['Re', title_str], 'fontsize',14);
grid on; xlim([w_min,w_max]);ylim([(min(Y_real)-0.5),(max(Y_real)+0.5)]);
xlabel(w_xLabelStr, 'fontsize',14);
line([w_min,w_max],[0,0], 'color','black');
set(gca,'xtick',w_x_tick);
set(gca,'xticklabel', w_x_tick_label_cell );

subplot(2,2,3); 
plot(w, Y_imag);title(['Im', title_str], 'fontsize',14);
grid on; xlim([w_min,w_max]);ylim([(min(Y_imag)-0.5),(max(Y_imag)+0.5)]);
xlabel(w_xLabelStr, 'fontsize',14);
line([w_min,w_max],[0,0], 'color','black');
set(gca,'xtick',w_x_tick);
set(gca,'xticklabel', w_x_tick_label_cell );


subplot(2,2,4); 
plot(w, Y_abs);title(['Mod', title_str], 'fontsize',14);
grid on; xlim([w_min,w_max]);ylim([(min(Y_abs)-0.5),(max(Y_abs)+0.5)]);
xlabel(w_xLabelStr, 'fontsize',14);
line([w_min,w_max],[0,0], 'color','black');
set(gca,'xtick',w_x_tick);
set(gca,'xticklabel', w_x_tick_label_cell );

function tick_str = get_tick_str(num_val, denorm_val, unit_str)

num_val    = round(num_val*1000); 
denorm_val = round(denorm_val*1000); 
gcd_n_d = gcd(num_val, denorm_val);
num_val = num_val/gcd_n_d;
denorm_val = denorm_val/gcd_n_d;

if(num_val == 0) 
  tick_str = ['0'];
else
  if (denorm_val == 1)
    denorm_str = unit_str;
  else
    denorm_str = [unit_str, '/', num2str(denorm_val)];
  end

  if(num_val == 1)
    num_str = [];
  elseif (num_val == -1)
    num_str = ['-'];
  else
    num_str = num2str(num_val);
  end
  tick_str = [num_str, denorm_str];
end



function ver_num = get_matlab_version()
ver_str_all = version; % such as 7.5.0.342 (R2007b)
ver_str = [ver_str_all(1), ver_str_all(2),ver_str_all(3)]; % such as 7.5
ver_num = str2num(ver_str);

function w_x_tick = get_w_x_tick_vec(w_min, w_max, w_x_tick_delta)

if(w_min < 0)
  w_x_tick_R = [0:w_x_tick_delta: w_max];
  temp       = [w_x_tick_delta:w_x_tick_delta:-w_min];
  w_x_tick_L = fliplr(-temp);
  w_x_tick   = [w_x_tick_L, w_x_tick_R];
else
  w_x_tick = [w_min:w_x_tick_delta:w_max];
end

if(w_x_tick(1) > w_min)
  w_x_tick = [w_min, w_x_tick];
end
if(w_x_tick(length(w_x_tick)) < w_max)
  w_x_tick = [w_x_tick, w_max];
end



function w_x_tick_label_cell = get_w_x_tick_label_cell(w_x_tick, w_x_tick_delta)
N_w_x_tick = length(w_x_tick);
w_x_tick_label_cell = cell(N_w_x_tick,1);
ver_num = get_matlab_version();
if(ver_num >= 8.4)
  tick_unit_str = '\pi';
else
  tick_unit_str = 'pi';
end
for(idx = 1:N_w_x_tick)
  num_val    = w_x_tick(idx)/w_x_tick_delta   ;
  denorm_val = pi/w_x_tick_delta   ;
  tick_str = get_tick_str(num_val, denorm_val, tick_unit_str);
  w_x_tick_label_cell{idx} = tick_str; 
end

2.3.2 实验结果

【FPGA作业】第二章
图2. 6 [-3,3], N=7
【FPGA作业】第二章
图2. 7 [0,6], N=7
【FPGA作业】第二章
图2. 8 [0,7], N=8
【FPGA作业】第二章
图2. 9 [-3,4], N=8
【FPGA作业】第二章
图2. 10 [-3,0] ,N=4
【FPGA作业】第二章
图2. 11 [-2,1], N=4
【FPGA作业】第二章
图2. 12 [-1,2], N=4
【FPGA作业】第二章
图2. 13 [0,3], N=4

2.3.3 实验分析

  • 对N点的矩形窗序列
  • 时域相位不同（即序列最左边和最右边的位置），DTFT模值大小不变。如[-3,3]和[0,6]的模值峰值均为7.4，主瓣、旁瓣宽度和旁瓣个数均相同。
  • 时域相位相同的DTFT实部主瓣比时域相位不同的实部主瓣窄，如[-3,3]和[0,6]，前者主瓣宽pi/2，后者pi/4；时域相位相同的DTFT虚部大小随频率增大减小，而时域相位不同的DTFT虚部大小为周期函数，如[0,7]和[-3,4]。
  • 偶对称序列，DTFT的虚部为0，如[-3,3]
  • 将序列按照纵轴翻转，DTFT的实部不变，虚部按照横轴翻转。

2.4 有限长正弦序列的DTFT

有限长序列相当于无限长序列与有限长矩形窗相乘得到，时域相乘，频域卷积

2.4.1 matlab程序

主程序，包括参数设置sin_dtft_plot.m

clc 
clear
close all
fig_fname = 'plot_dtft.png';

w_min = -1*pi;
w_max =  1*pi;

w_delta = pi/1000;

n_min =    0;
n_max =    15;

f0 = 125;
fs = 1000;

n = [n_min:1:n_max];
w = [w_min:w_delta:w_max];

N = length(n);

x = cos(  2*pi*f0/fs*n);

Y = func_calc_dtft(x, n, w); 
x_dscp_str = ['x[n] = sin(2\pi\cdotf0/fs\cdotn), f0=', num2str(f0),',fs=' num2str(fs)];
h = func_plot_dtft(x,n,w,Y, x_dscp_str);
dtft计算func_calc_dtft.m
% Y(w) = DTFT(x[n]);
function Y = func_calc_dtft(x, n, w) 

N     = length(n);
N_w   = length(w);

Y = zeros(N_w,1);
for idx = 1:N_w
  e_jw = exp(-j*w(idx)*n);
  Y(idx) = sum(e_jw .* x);
end

图像显示，包括信号时域图、dtft实部、dtft虚部和dtft模值func_plot_dtft.m

% file: func_plot_dtft.m
% test with Matlab 2007, 2014
% Y[w] = DTFT(x[n]);
% plot time sequence x[n], real/imag part of Y[w], abs of Y[w]
function h = func_plot_dtft(x,n,w,Y, x_dscp_str);
h = figure; 

N = length(n);
n_min = min(n);
n_max = max(n);
w_min = min(w);
w_max = max(w);
N_w = length(w);
w_delta = (w_max-w_min)/(N_w-1);

w_x_tick_delta = 1*pi/4 ;

w_x_tick = get_w_x_tick_vec(w_min, w_max, w_x_tick_delta);

w_x_tick_label_cell = get_w_x_tick_label_cell(w_x_tick, w_x_tick_delta);

ver_num = get_matlab_version();
if(ver_num >= 8.4)
  if(length(w_x_tick_label_cell) >= 12)
    x_tick_rot_angle = 270;
  else
    x_tick_rot_angle = 0;
  end
end

Y_real = real(Y); Y_imag = imag(Y); Y_abs = abs(Y);

subplot(2,2,1);
title_str = {x_dscp_str; ['n \in [',int2str(n_min),',', int2str(n_max), ... 
            '], ', int2str(N), ' samples']};
stem(n, x, 'fill'); title(title_str, 'fontsize',14);
y_lim_max = 1.1*max(abs(x)); y_lim_min = -y_lim_max; 
grid on; xlim([n_min-1,n_max+1]); ylim([y_lim_min, y_lim_max]);

subplot(2,2,2);
Left_str  = get_tick_str((w_min/pi),1,'\pi');
Right_str = get_tick_str((w_max/pi),1,'\pi');
w_xLabelStr =['\omega: [', Left_str,',', Right_str, ']'];
title_str = ['\{X(e^{j\omega})\}, \Delta\omega = \pi/', num2str(pi/w_delta)];
plot(w, Y_real);title(['Re', title_str], 'fontsize',14);
grid on; xlim([w_min,w_max]);ylim([(min(Y_real)-0.5),(max(Y_real)+0.5)]);
xlabel(w_xLabelStr, 'fontsize',14);
line([w_min,w_max],[0,0], 'color','black');
set(gca,'xtick',w_x_tick);
set(gca,'xticklabel', w_x_tick_label_cell );
if(ver_num >= 8.4) set(gca,'XTickLabelRotation',x_tick_rot_angle); end  

subplot(2,2,3); 
plot(w, Y_imag);title(['Im', title_str], 'fontsize',14);
grid on; xlim([w_min,w_max]);ylim([(min(Y_imag)-0.5),(max(Y_imag)+0.5)]);
xlabel(w_xLabelStr, 'fontsize',14);
line([w_min,w_max],[0,0], 'color','black');
set(gca,'xtick',w_x_tick);
set(gca,'xticklabel', w_x_tick_label_cell );
if(ver_num >= 8.4) set(gca,'XTickLabelRotation',x_tick_rot_angle); end  


subplot(2,2,4); 
plot(w, Y_abs);title(['Mod', title_str], 'fontsize',14);
grid on; xlim([w_min,w_max]);ylim([(min(Y_abs)-0.5),(max(Y_abs)+0.5)]);
xlabel(w_xLabelStr, 'fontsize',14);
line([w_min,w_max],[0,0], 'color','black');
set(gca,'xtick',w_x_tick);
set(gca,'xticklabel', w_x_tick_label_cell );
if(ver_num >= 8.4) set(gca,'XTickLabelRotation',x_tick_rot_angle); end  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tick_str = get_tick_str(num_val, denorm_val, unit_str)

num_val    = round(num_val*1000); 
denorm_val = round(denorm_val*1000); 
gcd_n_d = gcd(num_val, denorm_val);
num_val = num_val/gcd_n_d;
denorm_val = denorm_val/gcd_n_d;

if(num_val == 0) 
  tick_str = ['0'];
else
  if (denorm_val == 1)
    denorm_str = unit_str;
  else
    denorm_str = [unit_str, '/', num2str(denorm_val)];
  end

  if(num_val == 1)
    num_str = [];
  elseif (num_val == -1)
    num_str = ['-'];
  else
    num_str = num2str(num_val);
  end
  if (num_val > 0)
    num_str = ['+',num_str];
  end
  tick_str = [num_str, denorm_str];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ver_num = get_matlab_version()
ver_str_all = version; % such as 7.5.0.342 (R2007b)
ver_str = [ver_str_all(1), ver_str_all(2),ver_str_all(3)]; % such as 7.5
ver_num = str2num(ver_str);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function w_x_tick = get_w_x_tick_vec(w_min, w_max, w_x_tick_delta)

if(w_min < 0)
  w_x_tick_R = [0:w_x_tick_delta: w_max];
  temp       = [w_x_tick_delta:w_x_tick_delta:-w_min];
  w_x_tick_L = fliplr(-temp);
  w_x_tick   = [w_x_tick_L, w_x_tick_R];
else
  w_x_tick = [w_min:w_x_tick_delta:w_max];
end

if(w_x_tick(1) > w_min)
  w_x_tick = [w_min, w_x_tick];
end
if(w_x_tick(length(w_x_tick)) < w_max)
  w_x_tick = [w_x_tick, w_max];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function w_x_tick_label_cell = get_w_x_tick_label_cell(w_x_tick, w_x_tick_delta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N_w_x_tick = length(w_x_tick);
w_x_tick_label_cell = cell(N_w_x_tick,1);
ver_num = get_matlab_version();
if(ver_num >= 8.4)
  tick_unit_str = '\pi';
else
  tick_unit_str = 'pi';
end
for(idx = 1:N_w_x_tick)
  num_val    = w_x_tick(idx)/w_x_tick_delta   ;
  denorm_val = pi/w_x_tick_delta   ;
  tick_str = get_tick_str(num_val, denorm_val, tick_unit_str);
  w_x_tick_label_cell{idx} = tick_str; 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2.4.2 实验结果

【FPGA作业】第二章
图2. 14
模值曲线在对应的矩形窗零值栅格频率位置上的点，仍为0.

【FPGA作业】第二章
图2. 15
【FPGA作业】第二章
图2. 16
【FPGA作业】第二章
图2. 17
【FPGA作业】第二章
图2. 18
【FPGA作业】第二章
图2. 19
【FPGA作业】第二章
图2. 20
【FPGA作业】第二章
图2. 21
【FPGA作业】第二章
图2. 22
【FPGA作业】第二章
图2. 23
【FPGA作业】第二章
图2. 24
【FPGA作业】第二章
图2. 25
【FPGA作业】第二章
图2. 26

2.4.3 实验分析

  • 对实数序列信号，DTFT的实部和模值为偶函数，DTFT的虚部为奇函数
  • 正弦序列的DTFT曲线峰值，不一定出现在信号归一化数字角频率的位置上，因为分辨率导致的误差
  • 单音正弦信号频率过低时，由于靠近零点，在DTFT曲线上将无法分辨出正负频率峰值
  • 使用相同长度、不同区间的窗截取同一正弦信号时，信号经过DTFT后，模值保持不变，实部和虚部的过零点不变，但实部和虚部的数值大小改变。

2.5 DTFT和DFT关系探寻

DFT是0~2pi区间上对DTFT的等间隔采样，采样间隔大小为2pi/N

2.5.1 matlab程序

主程序，包括参数设置sin_dtft_dft_plot.m

clc 
clear
close all
fig_fname = 'plot_dtft.png';
% to plot DFT with DTFT, must set w range as [0,2pi]
w_min =  0;
w_max =  2*pi;

w_delta = pi/1000;

n_min =    0;
n_max =    16-1;

f0 = 10;
fs = 100;

n = [n_min:1:n_max];
w = [w_min:w_delta:w_max];

N = length(n);

x = sin(2*pi*f0/fs*n);

Y     = func_calc_dtft(x, n, w); 
Y_DFT = fft(x);

x_dscp_str = ['x[n] = sin(2\pi\cdotf0/fs\cdotn), f0=', num2str(f0),',fs=' num2str(fs)];
h = func_plot_dtft_dft(x,n,w,Y,Y_DFT, fs, x_dscp_str);
计算信号的dtft值func_calc_dtft.m
% Y(w) = DTFT(x[n]);
function Y = func_calc_dtft(x, n, w) 

N     = length(n);
N_w   = length(w);

Y = zeros(N_w,1);
for idx = 1:N_w
  e_jw = exp(-j*w(idx)*n);
  Y(idx) = sum(e_jw .* x);
end

绘图，包括信号时域图、dtft和dft的实部、虚部和模值func_plot_dtft_dft.m

% file: func_plot_dtft_dft.m
% test with Matlab 2007, 2014
% Y[w] = DTFT(x[n]);
% plot time sequence x[n], real/imag part of Y[w], abs of Y[w]
function h = func_plot_dtft_dft(x,n,w,Y, Y_DFT, fs, x_dscp_str);
h = figure; 

N = length(n);
n_min = min(n);
n_max = max(n);
w_min = min(w);
w_max = max(w);
if(w_min ~= 0   ) disp('WARNING, to plot DFT with DTFT w_min must be 0'); end
if(w_max ~= 2*pi) disp('WARNING, to plot DFT with DTFT w_max must be 2*pi'); end
N_w = length(w);
w_delta = (w_max-w_min)/(N_w-1);

w_x_tick_delta = 2*pi/16 ;

w_x_tick = get_w_x_tick_vec(w_min, w_max, w_x_tick_delta);

w_x_tick_label_cell = get_w_x_tick_label_cell(w_x_tick, w_x_tick_delta);

ver_num = get_matlab_version();
if(ver_num >= 8.4)
  if(length(w_x_tick_label_cell) >= 12)
    x_tick_rot_angle = 270;
  else
    x_tick_rot_angle = 0;
  end
end

Y_real = real(Y); Y_imag = imag(Y); Y_abs = abs(Y);

N_DFT = length(Y_DFT);
w_DFT = n*2*pi/N;

Y_DFT_real = real(Y_DFT); Y_DFT_imag = imag(Y_DFT); Y_DFT_abs = abs(Y_DFT);

subplot(2,2,1);
title_str = {x_dscp_str; ['n \in [',int2str(n_min),',', int2str(n_max), ... 
            '], ', int2str(N), ' samples']};
stem(n, x, 'fill'); title(title_str, 'fontsize',14);
y_lim_max = 1.1*max(abs(x)); y_lim_min = -y_lim_max; 
grid on; xlim([n_min-1,n_max+1]); ylim([y_lim_min, y_lim_max]);

subplot(2,2,2);
Left_str  = get_tick_str((w_min/pi),1,'\pi');
Right_str = get_tick_str((w_max/pi),1,'\pi');
w_xLabelStr =['\omega: [', Left_str,',', Right_str, ']'];
title_str = ['\{X(e^{j\omega})\}, \Delta\omega = \pi/', num2str(pi/w_delta)];
plot(w, Y_real);title(['Re', title_str], 'fontsize',14);
grid on; xlim([w_min,w_max]);ylim([(min(Y_real)-0.5),(max(Y_real)+0.5)]);
xlabel(w_xLabelStr, 'fontsize',14);
line([w_min,w_max],[0,0], 'color','black');
set(gca,'xtick',w_x_tick);
set(gca,'xticklabel', w_x_tick_label_cell );
if(ver_num >= 8.4) set(gca,'XTickLabelRotation',x_tick_rot_angle); end  
hold; stem(w_DFT, Y_DFT_real, 'filled');


subplot(2,2,3); 
plot(w, Y_imag);title(['Im', title_str], 'fontsize',14);
grid on; xlim([w_min,w_max]);ylim([(min(Y_imag)-0.5),(max(Y_imag)+0.5)]);
xlabel(w_xLabelStr, 'fontsize',14);
line([w_min,w_max],[0,0], 'color','black');
set(gca,'xtick',w_x_tick);
set(gca,'xticklabel', w_x_tick_label_cell );
if(ver_num >= 8.4) set(gca,'XTickLabelRotation',x_tick_rot_angle); end  
hold; stem(w_DFT, Y_DFT_imag, 'filled');


subplot(2,2,4); 
plot(w, Y_abs);title(['Mod', title_str], 'fontsize',14);
grid on; xlim([w_min,w_max]);ylim([(min(Y_abs)-0.5),(max(Y_abs)+0.5)]);
xlabel(w_xLabelStr, 'fontsize',14);
line([w_min,w_max],[0,0], 'color','black');
set(gca,'xtick',w_x_tick);
set(gca,'xticklabel', w_x_tick_label_cell );
if(ver_num >= 8.4) set(gca,'XTickLabelRotation',x_tick_rot_angle); end  
hold; stem(w_DFT, Y_DFT_abs, 'filled');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tick_str = get_tick_str(num_val, denorm_val, unit_str)

num_val    = round(num_val*1000); 
denorm_val = round(denorm_val*1000); 
gcd_n_d = gcd(num_val, denorm_val);
num_val = num_val/gcd_n_d;
denorm_val = denorm_val/gcd_n_d;

if(num_val == 0) 
  tick_str = ['0'];
else
  if (denorm_val == 1)
    denorm_str = unit_str;
  else
    denorm_str = [unit_str, '/', num2str(denorm_val)];
  end

  if(num_val == 1)
    num_str = [];
  elseif (num_val == -1)
    num_str = ['-'];
  else
    num_str = num2str(num_val);
  end
  if (num_val > 0)
    num_str = ['+',num_str];
  end
  tick_str = [num_str, denorm_str];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ver_num = get_matlab_version()
ver_str_all = version; % such as 7.5.0.342 (R2007b)
ver_str = [ver_str_all(1), ver_str_all(2),ver_str_all(3)]; % such as 7.5
ver_num = str2num(ver_str);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function w_x_tick = get_w_x_tick_vec(w_min, w_max, w_x_tick_delta)

if(w_min < 0)
  w_x_tick_R = [0:w_x_tick_delta: w_max];
  temp       = [w_x_tick_delta:w_x_tick_delta:-w_min];
  w_x_tick_L = fliplr(-temp);
  w_x_tick   = [w_x_tick_L, w_x_tick_R];
else
  w_x_tick = [w_min:w_x_tick_delta:w_max];
end

if(w_x_tick(1) > w_min)
  w_x_tick = [w_min, w_x_tick];
end
if(w_x_tick(length(w_x_tick)) < w_max)
  w_x_tick = [w_x_tick, w_max];
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function w_x_tick_label_cell = get_w_x_tick_label_cell(w_x_tick, w_x_tick_delta)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N_w_x_tick = length(w_x_tick);
w_x_tick_label_cell = cell(N_w_x_tick,1);
ver_num = get_matlab_version();
if(ver_num >= 8.4)
  tick_unit_str = '\pi';
else
  tick_unit_str = 'pi';
end
for(idx = 1:N_w_x_tick)
  num_val    = w_x_tick(idx)/w_x_tick_delta   ;
  denorm_val = pi/w_x_tick_delta   ;
  tick_str = get_tick_str(num_val, denorm_val, tick_unit_str);
  w_x_tick_label_cell{idx} = tick_str; 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2.5.2 实验结果

【FPGA作业】第二章
图2. 27
DFT的模值峰值刚好落在DTFT模值的峰值上
【FPGA作业】第二章
图2. 28
【FPGA作业】第二章
图2. 29
上面两个信号的DFT的模值峰值均未落在DTFT模值的峰值上

2.5.3 实验分析

假设有一个单音频率为f0的正弦信号采样序列，采样条件满足奈奎斯特定律。使用fs的采样率，记录了时间长度T、记录个数N个样点
* 如果用序列的DFT结果估计正弦信号的频率，采样个数越多的情况下，频率估计越准确
* 保持序列的时间长度T不变，提高采样率，频率分辨率提高
* 保持序列的记录个数N不变，提高采样率，频率分辨率提高
* 保持序列的采样率fs不变，提高记录个数N，频率分辨率提高
* 对序列补零可以使频率序列更为细致，但不能提高序列的频率分辨率
* 在序列后面补零与不补零得到的频谱一致，但更为细致；在序列前面补零，得到频谱的幅值相同，相位改变，相当于多采样了几个点

上一篇： hadoop小文件处理

下一篇： vector详解

推荐阅读