mydftdemo_1D.m

%==========================================================================
%                            一维信号DFT测试程序
%    dt:   时间采样间隔
%    L:    信号x时长
%    n:    信号x点数 (n = L/dt+1）
%    fs:   采样频率 (fs = 1/dt)
%    fn:   Nyquist频率（fn = fs/2=0.5*(1/dt)）
%    df:   频率精度(df = fs/(n-1)) 
%
%                 作者：彭真明
%                 单位：电子科技大学光电信息学院
%                 时间：2015.03.21
%==========================================================================
clc,clf,clear all, close all
%==========================================================================
% 初始参数设定
%==========================================================================
% 设定连续信号表达式为:
%      signal: x = sin(2*pi*10*t)+3*sin(2*pi*50*t)+2*sin(2*pi*120*t)
%==========================================================================
L = 0.6;       % 信号长度(秒)
Tmax = 1/10;   % 最大周期(秒)
fmax = 120;    % 单位：Hz
fn   = 2*fmax; % Nyquist频率
fsc  = fn;     % The critical sampling frequency.
dtc  = 1/fsc;  % The critical sampling time.
%==========================================================================
% 采样参数设定
%==========================================================================
dt = 0.0001, fs = 1./dt;
%==========================================================================
if fs < fsc
    warning('Smapling time is too small!!')
    return;
end
%==========================================================================



% 时间点赋值（单个周期）
t = 0:dt:Tmax;

% 按给定函数计算各时间点对应的x（即离散化处理）
x = sin(2*pi*10*t)+3*sin(2*pi*50*t)+2*sin(2*pi*120*t);

% 绘制y曲线（一个最大周期内）
subplot(311)
plot(1000*t,x)
title('Initial signal')
xlabel('time (ms)')

% 加随机噪声
y = x + 1.2*randn(size(t));

% 绘制y曲线（一个最大周期内）
subplot(312)
plot(1000*t,y)
title('Signal corrupted with zero-mean random noise')
xlabel('time (ms)')

N  = length(y);
df = fs/(N-1); % 频率间隔/精度

% DFT计算
%==========================================================================
% 变换前（时域）处理，实现频域的频谱中心化
for i=1:N
    y(i)=y(i)*(-1)^i;
end
%==========================================================================
y_dft = myDFT(y);
%f = df*(0:N-1);           % 对称频谱
%f = df*(0:N/2-1);         % 实频
f = df*(-(N-1)/2:(N-1)/2); % 对称频谱(含负频)
subplot(313)
plot(f,abs(y_dft));
title('Frequency spectrum of signal')
xlabel('frequency (Hz)')
%xlim([0 fs/2])

% set(gca,'XAxisLocation','top');   x坐标刻度在上
% set(gca,'YAxisLocation','right'); y坐标刻度在右