clear;clc;
Fs=120;                               % 采样频率 120 Hz, 大于两倍的最高频率
T=1/Fs;                               % 采样间隔
N=1024;                               % 采样点数
t=(0:1:N-1)*T;                        % 时间序列
t=t';                                 % 时间序列转置
 
% time-domain signal
y = 5 ...                             % 直流分量
    + (2)*sin(2*pi*20*t) ...          % 20 Hz
    + (3)*sin(2*pi*30*(t)) ...        % 30 Hz
    + (4)*sin(2*pi*40*(t));           % 40 Hz

%% 绘制时域信号
figure
plot(t,y)
xlabel('时间')
ylabel('信号值')
title('时域信号')

%% fft变换
Y=fft(y);          % Y为fft变换结果，复数向量
Y=Y(1:N/2+1);      % 只看变换结果的一半即可
A=abs(Y);          % 复数的幅值（模）
f=(0:1:N/2)*Fs/N;  % 生成频率范围
f=f';              % 转置成列向量
 
%% 幅值修正
A_adj=zeros(N/2+1,1);
A_adj(1)=A(1)/N;        % 频率为0的位置除以N
A_adj(end)=A(end)/N;    % 频率为Fs/2的位置也除以N
A_adj(2:end-1)=2*A(2:end-1)/N;  %其余频率点乘以2/N
 
%% 绘制频率幅值图
figure
subplot(2,1,1)
plot(f,A_adj)
xlabel('频率 (Hz)')
ylabel('幅值 (修正后)')
title('FFT变换幅值图')
grid on
 
%% 绘制频谱相位图
subplot(2,1,2)
phase_angle=angle(Y);   % angle函数的返回结果为弧度
phase_angle=rad2deg(phase_angle);%将弧度转化为度
plot(f,phase_angle)
xlabel('频率 (Hz)')
ylabel('相位角 (degree)')
title('FFT变换相位图')
grid on

import numpy as np
from scipy.fftpack import fft,fftshift
import matplotlib.pyplot as plt


N = 1024                        # 采样点数
sample_freq=120                 # 采样频率 120 Hz, 大于两倍的最高频率
sample_interval=1/sample_freq   # 采样间隔
signal_len=N*sample_interval    # 信号长度
t=np.arange(0,signal_len,sample_interval)

signal = 5 + 2 * np.sin(2 * np.pi * 20 * t) + 3 * np.sin(2 * np.pi * 30 * t) + 4 * np.sin(2 * np.pi * 40 * t)  # 采集的信号

fft_data = fft(signal)

# 在python的计算方式中，fft结果的直接取模和真实信号的幅值不一样。
# 对于非直流量的频率，直接取模幅值会扩大N/2倍， 所以需要除了N乘以2。
# 对于直流量的频率(0Hz)，直接取模幅值会扩大N倍，所以需要除了N。
fft_amp0 = np.array(np.abs(fft_data)/N*2)   # 用于计算双边谱
fft_amp0[0]=0.5*fft_amp0[0]
N_2 = int(N/2)
fft_amp1 = fft_amp0[0:N_2]  # 单边谱
fft_amp0_shift = fftshift(fft_amp0)    # 使用fftshift将信号的零频移动到中间

# 计算频谱的频率轴
list0 = np.array(range(0, N))
list1 = np.array(range(0, int(N/2)))
list0_shift = np.array(range(0, N))
freq0 = sample_freq*list0/N        # 双边谱的频率轴
freq1 = sample_freq*list1/N        # 单边谱的频率轴
freq0_shift=sample_freq*list0_shift/N-sample_freq/2  # 零频移动后的频率轴

# 绘制结果
plt.figure()
# 原信号
plt.subplot(221)
plt.plot(t, signal)
plt.title(' Original signal')
plt.xlabel('t (s)')
plt.ylabel(' Amplitude ')
# 双边谱
plt.subplot(222)
plt.plot(freq0, fft_amp0)
plt.title(' spectrum two-sided')
plt.ylim(0, 6)
plt.xlabel('frequency  (Hz)')
plt.ylabel(' Amplitude ')
# 单边谱
plt.subplot(223)
plt.plot(freq1, fft_amp1)
plt.title(' spectrum single-sided')
plt.ylim(0, 6)
plt.xlabel('frequency  (Hz)')
plt.ylabel(' Amplitude ')
# 移动零频后的双边谱
plt.subplot(224)
plt.plot(freq0_shift, fft_amp0_shift)
plt.title(' spectrum two-sided shifted')
plt.xlabel('frequency  (Hz)')
plt.ylabel(' Amplitude ')
plt.ylim(0, 6)

plt.show()