02-混合信號的頻譜分析

本文轉載自查看原文 2021-09-03 18:53 120

本文主要是分析由多個正弦或者余弦信號組成的混合信號，可否求解出各組成信號部分的幅值、角頻率及初始相位角，並進行驗證。

參考鏈接：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('初始相位值（角度）');