信号波形的产生与频谱的MATLAB程序.docx

目录
目录 ………………………………………………………………………………………………
正文…………………………………………………………………………………………………
1 基本波形………………………………………………………………………………
1.1 正弦波…………………………………………………………………………
1.2 三角波和锯齿波………………………………………………………………
1.3 方波和不对称方波……………………………………………………
1.4伪码……………………………………………………………………………………
2 调幅………………………………………………………………………………
2.1常规双边带调幅（AM）……………………………………………………………
2.2抑制载波双边带调幅(DSB)……………………………………………………………
2.3 单边带调幅（SSB）…………………………………………………………………
3.4振幅键控（ASK）……………………………………………………………………
3 调频………………………………………………………………………………
3.1单音频调频……………………………………………………………………………
3.2锯齿波调频……………………………………………………………………………
3.3 三角波调频…………………………………………………………………………
3.4 移频键控（FSK）………… ……………………………………………………………
4 调相……………………………………………………………………………………………
4.1 方波调相…………………………………………………………………………………
4.2 伪码调相…………………………………………………………………………………

一、 基本波形
（1）正弦波

为正弦波频率，为初始相位
理论频谱：
为正弦波角频率，且
时域波形：

图1-1-1 正弦波的时域波形（横坐标为时间t，纵坐标为幅度）
图1-1-1为正弦波的时域波形。从图中可以看出，正弦波的频率为50Hz，峰值为1，为1，所以初始相位为。
设采样率=1000Hz，一个周期采样20个点。

频域波形：

图1-1-2 正弦波的频域波形（横坐标为频率HZ，纵坐标为归一化幅度）
从图中可以看出，正弦波在正负50HZ的地方有二根谱线，与理论相符。
MATLAB程序：
% 正弦波波形
fm=50; %正弦波频率
fs=1000; %采样频率
dt=1/fs;
phase=0; %初相
t=0:dt:0.1;
figure(1)
y=sin(2*pi*fm*t+phase); %正弦波
plot(t,y);
grid on
xlabel('t');
ylabel('正弦波');
N=1024;
%Plot amplitude spectrum
sf=fft(y,N);
f=-fs/2:fs/N:(fs/2-fs/N);
figure(2)
plot(f,fftshift(abs(sf))/max(abs(sf)));
axis([-100 100 0 1]);
grid on;
title('Amplitude Spectrum of正弦波')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

（2）三角波和锯齿波（sawtooth）
三角波和锯齿波时域表达式：f(t)=2*fm*twidth-1 ,0≤t≤T2-2*fm*twidth+1+width1-width ,T2≤t≤T

为三角波频率。
周期三角波的傅立叶级数ft=4Eπ21∞1n2sin2nπ2cosnω1t
周期锯齿波的傅立叶级数ft=Eπ1∞(-1)n+11nsinnω1t
产生一个周期为T、幅度在-1到+1之间的周期性三角波信号。其中 表示最大幅度出现的位置：即在一个周期内，信号从t=0到×T时函数值从-1到+1线性增加，而从×T到T又是从+1到-1线性下降。取值在0 ~ 1之间。 

图1-2-1 三角波和锯齿波的时域波形（横坐标为时间t，纵坐标为幅度）
图1-2-1为三角波锯齿波的时域波形。从图中可以看出，频率为50Hz，
设采样率=10000Hz，一个周期采样200个点。
频域波形：

图1-2-2 三角波和锯齿波的时域波形（横坐标为频率HZ，纵坐标为归一化幅度）。
由上图可以看出三角波的频谱含有奇次谐波（50HZ的奇数倍）的余弦分量，幅度以 1n2的规律收敛。锯齿波的频谱只包含正弦分量，幅度以1n的规律收敛。
MATLAB程序：
% 锯齿波和三角波波形
fm=50; %锯齿波频率
fs=10000; %采样频率
dt=1/fs;
t=0:dt:0.1;
y=sawtooth(2*pi*fm*t,0.5); %三角波
z=sawtooth(2*pi*fm*t,0); %锯齿波
x=sawtooth(2*pi*fm*t,1); %锯齿波
figure(1);
subplot(3,1,1);
plot(t,y);
grid on;
xlabel('t');
ylabel('三角波');
subplot(3,1,2);
plot(t,z);
grid on;
xlabel('t');
ylabel('锯齿波1');
subplot(3,1,3);
plot(t,x);
grid on;
xlabel('t');
ylabel('锯齿波2');

N=2^20;
%Plot amplitude spectrum
Y=fft(y,N);
Z=fft(z,N);
X=fft(x,N);
f=-fs/2:fs/N:(fs/2-fs/N);
figure(2)
subplot(3,1,1);
plot(f,fftshift(abs(Y))/max(abs(Y)));
axis([-300 300 0 1]);
grid on;
title('Amplitude Spectrum of 三角波')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')
subplot(3,1,2);
plot(f,fftshift(abs(Z))/max(abs(Z)));
axis([-300 300 0 1]);
grid on;
title('Amplitude Spectrum of 锯齿波1')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')
subplot(3,1,3);
plot(f,fftshift(abs(X))/max(abs(X)));
axis([-300 300 0 1]);
grid on;
title('Amplitude Spectrum of 锯齿波2')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

（3）周期对称方波和不对称方波
时域表达式：(t)= 1,0≤t≤T*duty100 -1,T*duty100≤t≤T

频谱：
ft=2*duty100-1+4π1∞1nsin(nπduty100)cosnω1t
当
duty=50时，ft=4π1∞1nsin(nπ2)cosnω1t
duty=25时，ft=12+4π1∞1nsin(nπ4)cosnω1t

产生一个频率为、幅度为±1的周期性方波信号。其中表示占空比，即在信号的一个周期中正值所占的百分比。 

图1-3-1 周期方波的时域波形（横坐标为时间t，纵坐标为幅度）
由上图可以看出周期方波的频率为50HZ，占空比为50%，不对称周期方波的频率为50HZ，占空比为25%。

频域波形:

图1-3-2 周期方波的频域波形（横坐标为HZ，纵坐标为归一化幅度）
由图可以看出对称方波无直流分量，只含有奇次谐波的余弦分量，幅度以1n的规律收敛，不对称方波含有直流，幅度包络为抽样函数:4*duty100sa(nπduty100)。
采样频率为10000HZ，每个周期采样200个点。
MATLAB程序：
%方波
fm=50; %频率
fs=10000; %采样频率
dt=1/fs;
t=0:dt:0.1;
y=square(2*pi*fm*t,50);
z=square(2*pi*fm*t,25);
figure(1);
subplot(2,1,1);
plot(t,y);
axis([0 0.1 -1.5 1.5]);
xlabel('t');
ylabel('周期方波');
subplot(2,1,2);
plot(t,z);
axis([0 0.1 -1.5 1.5]);
xlabel('t');
ylabel('不对称周期方波');

N=2^20;
%Plot amplitude spectrum
Y=fft(y,N);
Z=fft(z,N);
f=-fs/2:fs/N:(fs/2-fs/N);
figure(2)
subplot(2,1,1);
plot(f,fftshift(abs(Y))/max(abs(Y)));
axis([-300 300 0 1]);
grid on;
title('Amplitude Spectrum of 周期方波')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')
subplot(2,1,2);
plot(f,fftshift(abs(Z))/max(abs(Z)));
axis([-300 300 0 1]);
grid on;
title('Amplitude Spectrum of 不对称周期方波')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

（4）伪码
伪随机序列是一种可重复产生的类似于噪声的序列，常用的伪随机序列包括m序列，Gold序列
m序列的产生（最长线性反馈移位寄存器序列）：
一般来说一个n级反馈移位寄存器可能产生的最长的周期为2n-1，反馈电路需要满足一定的条件才能得到最长序列。
m序列的性质：
1. 均衡性
在m序列的一个周期中，“0”“1”数目基本相等，“1”比“0”多一个。
2. 游程分布
序列中取值相等的那些相继的元素合称为一个“游程”，游程中元素的个数称为游程的长度。长度为k的游程占总游程数的2-k，且连0连1的游程数各占一半。

3. 相位相加特性
一个m序列Mp与其经任意延迟移位产生的另一个不同序列Mr模2相加，得到的仍是Mp的某次延迟移位序列Ms，即Mp⊕Mr=Ms。
4. 自相关函数
定义序列X=(x1，x2，… ，xn)，xi∈+1,-1的循环自相关为：
Rxk=1ni=1Nxix[i+k]N，若xi=1-2ai，其中ai是m序列的输出，则xixj=1-2(xi⊕xj)，由于m序列的移位相加仍得到m序列，且m序列中1的个数比0的个数多1，应此
Rxk=1 , k=0-1N ,k=1,2…N-1
当N很大时，m序列的自相关趋于冲激函数，加上脉冲成形后，m序列信号是一个周期为Tc的周期函数
mt=1nj=-∝∝i=1Nxip[t-i-1Tc-jNTc]
其自相关函数也为周期函数
Rm(τ)=1NTc0NTcmtmt+τdt
求出：
R(τ)=1-N-1NTcτ-iNTc 0≤τ-iNTc≤Tc-1N others

对上式做傅立叶变换，求出其功率谱密度：
Ps(ω)=N+1N2[sinωTc2ωTc2]2n=-∞n≠0∞δω-2nπNTc+1N2δ(ω)
用MATLAB产生m序列，其生产多项式为：
gx=1+x2+x3

图1-4-1 伪随机序列的时域波形（横坐标为抽样点数，纵坐标为信号幅度）
采样频率为15000HZ，信号频率为150hz，每个周期采样100个点。由上图可以看出每个周期随机序列为000100110101111(m序列的0对应1，1对应-1)，上图显示了二个周期。

图1-4-2 伪随机序列的频域波形（横坐标为HZ，纵坐标幅度/db）
可以看出第一个零点位置是150HZ。

图1-4-3 伪随机序列的频域波形的局部放大图(横坐标为HZ，纵坐标幅度）可以看出信号在2nπNTc，即nNf(-∞≤n≤∞，f=150NHZ) 处有谱线（N=15），幅度满足功率谱密度表达式。

MATLAB程序：
%m序列发生器
clear all;
close all;
fm=150; %锯齿波频率
fs=15000; %采样频率
NUM=fs/fm;

g=19; %10011
state=8; %1000
DataLen=200;

%m序列产生
N=15;
mq=mgen(g,state,DataLen);

x_t=upsample(mq,NUM);
for k=1:NUM:NUM*DataLen-1
for i=1:NUM-1
x_t(k+i)=x_t(k);
end
end
x_t=2*x_t-1;
n=1:length(x_t);
figure(1)
plot(n,x_t);
grid on;
axis([0 3000 -2 2]);

N=2^20;
%Plot amplitude spectrum
sf=fft(x_t,N);
f=-fs/2:fs/N:(fs/2-fs/N);
figure(2)
plot(f,fftshift(abs(sf))/max(abs(sf)));
%plot(f,10*log10(fftshift(abs(sf))));
grid on;
axis([-600 600 0 1]);
%axis([-300 300 -30 50]);
title('Amplitude Spectrum of m序列')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

二、 调幅
（1）常规双边带调幅（AM）
设均值为零的模拟基带信号为m(t)，双边带调幅（AM）信号为：
sAMt=[A+mt]cos⁡(2*π*fc*t)
其中A是一常数，当A>mt时，称此调幅信号欠调幅；当A<mt时，称过调幅。当mt的频宽远小于载波信号时，欠调幅信号可以用包络检波的方式解调，而过调幅只能通过相关解调。
如果m(t)是确知信号，则频谱为：
Sf=12Mf-fc+Mf+fc+A2[δf-fc+δ(f+fc)] 设A=1
用Matlab产生一个频率为100kHz的余弦信号m(t)，并设fc=3MHZ，采样频率为30MHZ。
时域波形：

图2-1-1 AM调制信号时域波形（横坐标为时间t，纵坐标为幅度）
图中A=1，mt的幅度也为1。

图2-1-2 AM调制信号频域波形（横坐标为频率f，纵坐标为幅度）

图2-1-3 AM调制信号频域波形的放大（横坐标为频率f，纵坐标为幅度）
由频谱图可以看出载波频率fc=3MHZ，余弦调制信号频率为100kHz。

%显示AM调制波形
clear all;
A=1;
fs=30*10^6; %采样频率
dt=1/fs; %采样间隔
fc=3*10^6; %载波频率
fm=100*10^3; %信源频率
T=3*10^-5; %信号时长
t=0:dt:T;
mt=cos(2*pi*fm*t); %信源

%AM modulation
s_AM=(mt+A).*cos(2*pi*fc*t);
figure(1)
plot(t,s_AM);
hold on;
plot(t,mt+A,'r--');
title('AM调制信号');
xlabel('t');

N=2^10;
%Plot amplitude spectrum
sf=fft(s_AM,N);
f=-fs/2:fs/N:(fs/2-fs/N);
figure(2)
plot(f,fftshift(abs(sf))/max(abs(sf)));
grid on;
axis([-4*10^6 4*10^6 0 1]);
title('Amplitude Spectrum of AM调制信号')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')
figure(3)
plot(f,fftshift(abs(sf))/max(abs(sf)));
grid on;
axis([2*10^6 4*10^6 0 1]);
title('Amplitude Spectrum of AM调制信号')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')
（2）抑制载波双边带调幅(DSB)
设均值为零的模拟基带信号为m(t)，抑制载波双边带调幅(DSB)信号为：
sDSBt=mtcos⁡(2*π*fc*t)
如果m(t)是确知信号，则频谱为：
Sf=12[Mf-fc+Mf+fc]，其中Mf是m(t)的频谱。
相干解调原理：
rt=st*cos2*π*fc*t=12mt+12mt*cos⁡(4*π*fc*t)
再用低通滤波器将高频分量滤掉，就可恢复原始信息。
用Matlab产生一个频率为100kHz的余弦信号m(t)，并设fc=3MHZ，采样频率为30MHZ。
时域调制波形：

图2-2-1 DSB调制信号时域波形（横坐标为时间t，纵坐标为幅度）

图2-2-2 DSB调制信号频域波形（横坐标为频率f，纵坐标为幅度）

图2-2-3 DSB调制信号频域波形的放大（横坐标为频率f，纵坐标为幅度）
由频谱图可以看出载波频率fc=3MHZ，但是被抑制掉了，余弦调制信号频率为100kHz。
%显示DSB调制波形
clear all;
fs=30*10^6; %采样频率
dt=1/fs; %采样间隔
fc=3*10^6; %载波频率
fm=100*10^3; %信源频率
T=2*10^-5; %信号时长
t=0:dt:T;
mt=cos(2*pi*fm*t); %信源
% DSB modulation
s_DSB=mt.*cos(2*pi*fc*t);
figure(1)
plot(t,s_DSB);
hold on;
plot(t,mt,'r--');
title('DSB调制信号');
xlabel('t');

N=2^20;
%Plot amplitude spectrum
sf=fft(s_DSB,N);
f=-fs/2:fs/N:(fs/2-fs/N);
figure(2)
plot(f,fftshift(abs(sf))/max(abs(sf)));
grid on;
axis([-4*10^6 4*10^6 0 1]);
title('Amplitude Spectrum of DSB调制信号')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

figure(3)
plot(f,fftshift(abs(sf))/max(abs(sf)));
grid on;
axis([2*10^6 4*10^6 0 1]);
title('Amplitude Spectrum of AM调制信号')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

（3）单边带调幅（SSB）
设均值为零的模拟基带信号为m(t)，抑制载波双边带调幅(SSB)上边带信
号为：
SSSBt=mtcos2*π*fc*t-mtsin2*π*fc*t
=Re[(mt+jmt)ej2πfct]
=12[mt+jmtej2πfct+(mt-jmt)e-j2πfct]
故：Sf=12[M+f-fc+M-f+fc]
其中，m(t)是mt的希尔伯特变换。
如果 M(ω)是 m(t)的频谱，W(ω)是m(t)的频谱，有：
W(ω)=Mω[-jsgn(ω)]
用Matlab产生一个频率为100kHz的余弦信号m(t)，并设fc=3MHZ，采样频率为30MHZ。
时域波形：

图2-3-1 SSB调制信号时域波形（横坐标为时间t，纵坐标为幅度）

图2-3-2 SSB调制信号频域波形（横坐标为频率f，纵坐标为幅度）

图2-3-3 SSB调制信号频域波形的放大（横坐标为频率f，纵坐标为幅度）
由频谱图可以看出载波频率fc=3MHZ，但是被抑制掉了，余弦调制信号频率为100kHz，并且只有上边带。

MATLAB程序：
%显示SSB调制波形
clear all;
fs=30*10^6; %采样频率
dt=1/fs; %采样间隔
fc=3*10^6; %载波频率
fm=100*10^3; %信源频率
T=2*10^-5; %信号时长
t=0:dt:T;
mt=cos(2*pi*fm*t); %信源

%SSB modulation
s_SSB=real(hilbert(mt).*exp(j*2*pi*fc*t));
figure(1)
plot(t,s_SSB);
hold on;
plot(t,mt,'r--');
title('SSB调制信号');
xlabel('t');

N=2^20;
%Plot amplitude spectrum
sf=fft(s_SSB,N);
f=-fs/2:fs/N:(fs/2-fs/N);
figure(2)
plot(f,fftshift(abs(sf))/max(abs(sf)));
grid on;
axis([-4*10^6 4*10^6 0 1]);
title('Amplitude Spectrum of SBB调制信号')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

（4）二进制振幅键控（2ASK）
如果将二进制码元“0”对应信号0，“0”对应信号Acos2*π*fc*t，则2ASK信号可以写成如下表达式：
S2ASKt=ang(t-nTs)Acos2*π*fc*t
其中an∈0，1，gt=1 0≤t≪Ts0 others 。

可以看出，上式是数字基带信号mt=ang(t-nTs)经过DSB调制后形成的信号，其功率谱密度：

Psf=A24[Pmf-fc+Pmf+fc]，不妨设A=1。

用MATLAB产生独立等概率的二进制信源，频率为300kHz的余弦信号m(t)，并设fc=3MHZ，采样频率为30MHZ。
画出信号波形及功率谱如下：

图2-4-1 ASK调制信号时域波形（横坐标为时间t，纵坐标为幅度）
二进制码元的产生是用的第一节中的m序列（见图1-4-1）。

图2-4-2 ASK调制信号频域波形（横坐标为频率f，纵坐标为幅度）

图2-4-3 ASK调制信号频域波形的局部放大图（1）
从频域可以ASK调制就是将调制信号的基带频谱线性搬至中频处。第一零点的位置就是调制信号的频率300KHZ。

图2-4-3 ASK调制信号频域波形的局部放大图（2）
由上图看出在带宽内的离散谱线分布在f=300NKHZ（N=15），即20KHZ处。
虽然ASK调制是一种抑制载波的调制方式，但由于调制信号中的m序列没有进行矩形成形，其频谱中是含有直流分量的，导致调制后的载波分量还是比较大的。
程序：
%显示ASK调制波形
clear all;
fs=30*10^6; %采样频率
Ts=1/fs; %采样间隔
fc=3*10^6; %载波频率
fm=300*10^3; %信源频率
NUM=fs/fm;

g=19; %10011
state=8; %1000
DataLen=200;

%m序列产生
N=15;
mt=mgen(g,state,DataLen);

x_t=upsample(mt,NUM);
for k=1:NUM:NUM*DataLen-1
for i=1:NUM-1
x_t(k+i)=x_t(k);
end
end

n=1:length(x_t);
s_ASK=x_t.*cos(2*pi*fc/fs*n);
figure(1)
plot(n,s_ASK);
hold on;
plot(n,x_t,'r--');
axis([0 3000 -1 1]);

title('ASK调制信号');
xlabel('');

N=2^20;
%Plot amplitude spectrum
sf=fft(s_ASK,N);
f=-fs/2:fs/N:(fs/2-fs/N);
figure(2)
plot(f,fftshift(abs(sf))/max(abs(sf)));
%plot(f,10*log10(fftshift(abs(sf))));
grid on;
axis([-4*10^6 4*10^6 0 1]);
title('Amplitude Spectrum of ASK调制信号')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

三、 调频
当载波的频率变化与输入的基带信号幅度的变化成线性关系时，就构成了调频信号，调频信号可以写成：
st=Acos(2πfct+2πKf-∞tm(τ)dτ)，其中Kf称为调频偏移常数。
调频信号的频谱与输入信号的频谱之间不再是频率搬迁的关系，应此通常无法写出调频信号频谱的明确表达式，但调频信号的功率带宽与调频指数和输入信号的带宽有关，根据经验公式近似计算：
B=2(mf+1)fm，其中mf=∆fmaxfm，∆fmax和fm为最大频偏和信号带宽。

（1） 单音频调频
设mt=cos⁡(ωmt)，则st=Acos(2πfct+Kffmsin⁡(ωmt))
设mf=Kffm
st=Acos2πfct*cosmf*sinωmt
-Asin2πfct*sinmf*sinωmt
cosmf*sinωmt=Jomf+n=1∞2J2nmfcos⁡(2nωmt)
sinmf*sinωmt=n=1∞2J2n-1mfcos⁡((2n-1）ωmt)
其中Jnmf为第一类n阶贝塞尔函数。
频域表达式：
Sω=πAn=-∞∞Jnmf[δω-ωc-nωm+δ(ω+ωc+nωm)]

用Matlab产生一个频率为300kHz的余弦调制信号m(t)，并设fc=3MHZ，采样频率为30MHZ。mf为不同值时的波形：

图3-1-1 mf=0时，即无调制信号时的波形与频谱
无调制时Jomf=1，其他都为0，故只在fc处有谱线。

图3-1-2 mf=1的波形与频谱

图3-1-3 mf=1的频谱放大图，谱线分布在fc±nfm处，谱线幅度与
贝塞尔函数Jnmf有关。

图3-1-4 mf=5时的波形与频谱

图3-1-5 mf=5的频谱放大图，谱线分布在fc±nfm处，谱线幅度与
贝塞尔函数Jnmf有关。
由上图的比较可以看出，调频信号的一个特点是，它不改变总功率，而只是改变功率在各分量之间的分配关系。当调频偏移常数变大时，调频指数变大，调频信号带宽也变大，带宽为B=2(mf+1)fm。

%显示调频信号波形
clear all;
fs=30*10^6; %采样频率
dt=1/fs; %采样间隔
fc=3*10^6; %载波频率
fm=300*10^3; %信源频率
T=1*10^-5; %信号时长
kf=1500*10^3;
t=0:dt:T;
mt=cos(2*pi*fm*t); %调制信号

int_m(1)=0;
for i=1:length(t)-1
int_m(i+1)=int_m(i)+mt(i)*dt;
end

s_FM=cos(2*pi*fc*t+2*pi*kf*int_m);
figure(1)
subplot(2,1,1)
plot(t,s_FM);
hold on;
plot(t,mt,'r--');
title('单音调频信号');
xlabel('t');

N=2^20;
%Plot amplitude spectrum
sf=fft(s_FM,N);
f=-fs/2:fs/N:(fs/2-fs/N);
subplot(2,1,2)
plot(f,fftshift(abs(sf))/max(abs(sf)));
grid on;
axis([-10*10^6 10*10^6 0 1]);
title('Amplitude Spectrum of 单音调频信号')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

figure(2)
plot(f,fftshift(abs(sf))/max(abs(sf)));
grid on;
axis([-10*10^6 10*10^6 0 1]);
title('Amplitude Spectrum of 单音调频信号')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

（2）锯齿波调频
设mt为锯齿波，则st=Acos(2πfct+2πKf-∞tm(τ)dτ)
用Matlab产生一个频率为300kHz的余弦调制信号m(t)，并设fc=3MHZ，采样频率为30MHZ。mf为不同值时的波形：

图3-2-1 mf=1时的波形与频谱

图3-2-2 mf=5时的波形与频谱的放大图

图3-2-3 mf=5时的波形与频谱

图3-2-4 mf=5时的波形与频谱的放大图
分析结果与单音调频类似。

MATLAB程序：
%显示调频信号波形
clear all;
fs=30*10^6; %采样频率
dt=1/fs; %采样间隔
fc=3*10^6; %载波频率
fm=300*10^3; %信源频率
T=1*10^-5; %信号时长
kf=300*10^3;
t=0:dt:T;
mt=sawtooth(2*pi*fm*t,1); %锯齿波

int_m(1)=0;
for i=1:length(t)-1
int_m(i+1)=int_m(i)+mt(i)*dt;
end

s_FM=cos(2*pi*fc*t+2*pi*kf*int_m);
figure(1)
subplot(2,1,1);
plot(t,s_FM);
hold on;
plot(t,mt,'r--');
title('锯齿波调频信号');
xlabel('t');

N=2^20;
%Plot amplitude spectrum
sf=fft(s_FM,N);
f=-fs/2:fs/N:(fs/2-fs/N);
subplot(2,1,2);
plot(f,fftshift(abs(sf))/max(abs(sf)));
grid on;
axis([-10*10^6 10*10^6 0 1]);
title('Amplitude Spectrum of 锯齿波调频信号')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

figure(2)
plot(f,fftshift(abs(sf))/max(abs(sf)));
grid on;
axis([-10*10^6 10*10^6 0 1]);
title('Amplitude Spectrum of 锯齿波调频信号')
xlabel('Frequency (Hz)')
ylabel('|F(f)|')

（3）三角波调频
mt为三角波，画出波形：

图3-3-1 mf=1时的波形与频谱

图3-3-2 mf=1时的波形与频谱的放大图

图3-3-3 mf=5时的波形与频谱

图3-3-4 mf=5时的波形与频谱的放大图
分析原理和结果与单音调频类似。

%显示调频信号波形
clear all;
fs=30*10^6; %采样频率
dt=1/fs; %采样间隔
fc=3*10^6; %载波频率
fm=300*10^3; %信源频率
T=1*10^-5; %信号时长
kf=300*10^3;
t=0:dt:T;
mt=sawtooth(2*pi*fm*t,0.5); %三角波

int_m(1)=0;
for i=1:length(t)-1
int_m(i+1)=int_m(i)+mt(i)*dt;
end

s_FM=cos(2*pi*fc*t+2*pi*kf*int_m);
figure(1)
subplot(2,1,1);
plot(t,s_FM);
hold on;
plot(t,mt,'r--');
title('三角波调频信号');
xlabel('t');

N=2^20;
%Plot amplitude spectrum
sf=fft(s_FM,N);
f=-fs/2:fs/N:(fs/2-fs/N);
subplot(2,1,2);
plot(f,fftshift(abs(sf))/max(abs(sf)));
grid on;
axis([-10*10^6 10*10^6 0 1]