02-混合信号的频谱分析

本文主要是分析由多个正弦或者余弦信号组成的混合信号，可否求解出各组成信号部分的幅值、角频率及初始相位角，并进行验证。

参考链接：https://zhuanlan.zhihu.com/p/401208432

　　　https://mp.weixin.qq.com/s/MjA2wDENii5XGdUwFUHHdg

例1：

混合信号f(t)由N个余弦信号组成，其各部分对应的参数见表1所示；

clear all; close all;

A1 =5;                    %频率F1信号的幅度
A2 =10;                 %频率F2信号的幅度
A3 =15;                    %频率F2信号的幅度
A4 =20;                    %频率F2信号的幅度

F1 = 100;                   %信号1频率(Hz)
F2 = 1000;                   %信号2频率(Hz)
F3 = 500;
F4 = 800;

P1 = 0;                   %信号1相位(度)
P2 = 45;                    %信号2相位(度)
P3 = 90;
P4 = 135;

Fs = 5120;                  %采样频率(Hz)

N = 256;                     %采样点数，请注意时域采样N点，FFT时也是N点！
t = [0 : 1/Fs : N/Fs];        %采样时刻，注意不是序列是真正的时间！t = [0:N]*Ts

% 生成信号
signal1=A1*cos(2*pi*F1*t+pi*P1/180); % 周期T1 = 2*pi/(2*pi*F1) = 0.01s，采样周期Ts = 0.00019531 < (T1)/2
signal2=A2*cos(2*pi*F2*t+pi*P2/180); % 周期T2 = 2*pi/(2*pi*F2) = 0.001s，采样周期Ts = 0.00019531 < (T2)/2
signal3=A3*cos(2*pi*F3*t+pi*P3/180); 
signal4=A4*cos(2*pi*F4*t+pi*P4/180); 

signal5=A3*sin(2*pi*F3*t+pi*P3/180); 
signal6=A4*sin(2*pi*F4*t+pi*P4/180); 

% 原始信号+随机噪音
Y1 = wgn(1,length(signal1),0) + awgn(signal1,10,0);
Y2 = wgn(1,length(signal2),2)+ awgn(signal2,4,5,'linear');
Y3 = wgn(1,length(signal3),3) + awgn(signal3,4,'measured','linear'); % 
Y4 = wgn(1,length(signal4),4) + signal4;

% 合成信号
S1 = signal1 + signal2 + signal3 + signal4;
% 加噪音的合成信号
% S = Y1 + Y2 + Y3 + Y4;
S = signal1 + signal2 + signal5 + signal6;


%显示原始信号
figure(1);
% plot(t,S1,'b', t, S, 'r');
plot(S);
title('混合余弦+正弦信号');
legend('原始的混合信号', '加噪音的信号');
xlabel('时间');
ylabel('幅值');

%% 频率与幅值
Y = fft(S, N);
Ayy = (abs(Y));

% 因为MATLAB中FFT的变换矩阵不是一个酉矩阵(Unitary Matrix)，该阵除以1/sqrt(N)就是个酉矩阵。故经过变换后对信号有放大作用，
% 所以要在fft处理后结果除以N/2来，修正此“放大”作用。但是结果在直流的那一点是错误的，实际上直流应该除以N修正。
F = ([1:N]-1)*Fs/N;     % 换成实际的频率值， N*Fs/2 
Ayy = Ayy / (N /2);      %换算成实际的幅度
Ayy(1) = Ayy(1) /2; 

figure(2);
plot(F(1:N/2),Ayy(1:N/2)); 
title('实际幅度-频率曲线图');
xlabel('频率f/Hz');
ylabel('幅度值');

%% 频率与相位
% 相位的计算可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角度值，范围从-pi到pi。
% atan2(500, 148)=x,结果是弧度，换算为角度就是180*(-x)/pi=相位值。
Pyy = [1 : N/2];
for i = 1 : N/2
    Pyy(i) =phase(Y(i));               %计算相位
    Pyy(i) = Pyy(i) * 180 /pi;          %换算为角度
end;
figure(3);
plot(F(1 : N/2), Pyy(1 : N/2));      %显示相位图
title('相位-频率曲线图');
xlabel('频率f/Hz');
ylabel('初始相位值（角度）');