AM信号的调制与解调

1
2
3
4
5
6
7
8
9
10
11
clear;clc;close all;
echo on                   % 打开所有命令文件的显示方式，从这往后每一步的运行都会显示
t0 = .15;                   % signal duration 
ts = 0.001;                 % sampling interval 
fc = 250;                   % carrier frequency   
snr = 10;                   % SNR in dB (logarithmic) 
a = 0.85;                   % modulation index--调制指数
fs = 1 / ts;                  % sampling frequency     
t = [0 : ts : t0];              % time vector    
df = 0.2;                   % required frequency resolution--所需的频率分辨率
snr_lin = 10^(snr / 10);      % SNR 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
% 调制信号(m)的产生，AM已调信号(u)的产生
% 调制信号m 是一个分段表示的信号，在(0 , t0/3)上为 1  ，在(t0/3 , 2*t0/3)上为 -2 ，其他为0 
m = [ones(1,t0 / (3 * ts)),-2 * ones(1,t0 / (3 * ts)),zeros(1,t0 / (3 * ts) + 1)];
c = cos(2 * pi * fc .* t);                  % carrier signal
m_n = m / max(abs(m));                      % normalized message signal 归一化处理的调制信号
u = (1 + a * m_n) .* c;                     % modulated signal
% A0是加的直流，不携带信息，在此直接略去
subplot(2,2,1)
plot(t,m(1:length(t)))
axis([0 0.15 -2.1 2.1])
xlabel('Time')
title('The message signal')
subplot(2,2,2)
plot(t,c(1:length(t)))
axis([0 0.15 -2.1 2.1])
xlabel('Time')
title('The carrier') 
subplot(2,2,3)
plot(t,u(1:length(t)))
axis([0 0.15 -2.1 2.1])
xlabel('Time')
title('The modulated signal')

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
% 调制信号和已调信号的频谱特征，信号的功率
[M,m,df1] = fftseq(m,ts,df);               % Fourier transform
M = M / fs;                             % scaling(缩放)
f = [0 : df1 : df1 * (length(m)-1)] - fs / 2;          % frequency vector
[U,u,df1] = fftseq(u,ts,df);                 % Fourier transform
U = U / fs;                              % scaling
signal_power = spower(u(1 : length(t)));       % power in modulated signal
pmn = spower(m(1 : length(t))) / (max(abs(m)))^2;
eta = (a^2 * pmn) / (1 + a^2 * pmn);
% a = (Am/A0)^2 ，A0是AM信号加的直流信号
% pmn = ( avg( m(t) ) / Am )^2      Am是
% 化简后 eta = ( avg( m(t) )^2 ) / [ avg( m(t) )^2 + A0^2 ] = η （调制效率）
figure(2);
subplot(2,1,1);
plot(f,abs(fftshift(M))); 
xlabel('Frequency');
title('Spectrum of the message signal');     % 调制信号的频谱
subplot(2,1,2);
plot(f,abs(fftshift(U))); 
title('Spectrum of the modulated signal');
xlabel('Frequency');

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
% 噪声及叠加噪声后的信号的特征
noise_power = eta * signal_power / snr_lin;    % noise power
noise_std = sqrt(noise_power);            % noise standard deviation噪声标准差
noise = noise_std * randn(1,length(u));       % generate noise
r = u + noise;                           % add noise to the modulated signal
[R,r,df1] = fftseq(r,ts,df);                 % Fourier transform 
R = R / fs;
figure(3);
subplot(2,1,1);
plot(t,noise(1:length(t))) ;
title('noise sample') ;
xlabel('Time');
subplot(2,1,2);
plot(t,r(1:length(t))) ;
title('Signal and noise');
xlabel('Time');
figure(4);
subplot(2,1,1);
plot(f,abs(fftshift(U))) ;
title('Signal spectrum');
xlabel('Frequency');
subplot(2,1,2);
plot(f,abs(fftshift(R))) ;
title('Signal and noise spectrum');
xlabel('Frequency');

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
% AM的包络检波
env_r = env_phas(r);                      % envelope, when noise is present 当噪声存在时的包络
dem2 = 2 * (env_r - 1) / a;                     % demodulate in the presence of noise 当噪声存在时的解调
figure(5);
subplot(3,1,1);
plot(t,m(1 : length(t)));
axis([0 0.15 -2.1 2.1]);
xlabel('Time');
title('The message signal');
subplot(3,1,2);
plot(t, env_r (1:length(t)));
xlabel('Time');
title('The envelope of the modulated signal');
subplot(3,1,3);
plot(t, dem2 (1:length(t)));
xlabel('Time');
title('The demodulated signal in the presence of noise');

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
function [M,m,df]=fftseq(m,ts,df) 
%       [M,m,df]=fftseq(m,ts,df)
%       [M,m,df]=fftseq(m,ts)
%FFTSEQ     generates M, the FFT of the sequence m.     
%       The sequence is zero padded to meet the required frequency resolution df.
%       零填充序列，以满足所需的（输入）频率分辨率df
%       ts is the sampling interval. The output df is the final frequency resolution.
%       ts是采样间隔， 输出的df是 最终的 频率分辨率
%       Output m is the zero padded version of input m. M is the FFT.
%       输出m 是 输入m 的零填充版本                        M是FFT变化
fs = 1 / ts;
% nargin 针对当前正在执行的函数，返回函数调用中给定函数输入参数的数目。
if nargin == 2          %  [M,m,df]=fftseq(m,ts)
  n1 = 0;
else                    %  [M,m,df]=fftseq(m,ts,df)
  n1 = fs / df;
end
n2 = length(m);
% y = nextpow2(x);   2^y 为 大于等于x 的最小的二的整数次幂的数字 nextpow2(4)=2  nextpow2(5)=3
n = 2^(max(nextpow2(n1),nextpow2(n2)));
M = fft(m,n);       % Y = fft(X,n)
% 如果 X 是向量且 X 的长度小于 n，则为 X 补上尾零以达到长度 n
% 如果 X 是向量且 X 的长度大于 n，则对 X 进行截断以达到长度 n
m = [m,zeros(1,n - n2)];
df = fs / n;

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function [v,phi]=env_phas(x,ts,f0)
%		[v,phi]=env_phas(x,ts,f0)
%		      v=env_phas(x,ts,f0)
% ENV_PHAS   	Returns the envelope and the phase of the bandpass signal x.
              % 返回带通信号的 包络和相位
%		f0 is the center frequency. 
%		ts is the sampling interval.

% nargin 针对当前正在执行的函数，返回函数调用中给定函数输出参数的数目。
if nargout == 2
  z = loweq(x,ts,f0);   % 返回 信号x 的低通等效值
  phi = angle(z);
end
v = abs(hilbert(x));    % 希尔伯特变换，位移 π/2

1
2
3
4
5
6
7
8
9
10
11
function xl=loweq(x,ts,f0)
%       xl=loweq(x,ts,f0)
%LOWEQ      returns the lowpass equivalent of the signal x
          % 返回 信号x 的低通等效值
%       f0 is the center frequency. 
%       ts is the sampling interval.
   
t = [0 : ts : ts * (length(x)-1)];  
z = hilbert(x);
xl = z.* exp(-j * 2 * pi * f0 * t);

1
2
3
4
5
6
7
8
function p=spower(x)
%       p=spower(x)
%SPOWER     returns the power in signal x
          % 返回 信号x 的功率
p = (norm(x)^2) / length(x);
% norm(x,1),x是一个向量，norm是对向量中所有值的绝对值求和
% p就是x的平均功率