《信号与系统》- 第四章 傅里叶变换和系统的频域分析

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 生成方波信号
def square_wave(t, T):
    """生成周期为T的方波信号"""
    return np.where(np.mod(t, T) < T/2, 1, -1)

# 方波的傅里叶级数展开
def fourier_series_square(t, T, n_terms):
    """使用n_terms项傅里叶级数逼近方波"""
    result = 0
    for n in range(1, n_terms+1, 2):  # 只包含奇次谐波
        result += (4/(n*np.pi)) * np.sin(2*np.pi*n*t/T)
    return result

# 时间范围
t = np.linspace(-2, 2, 1000)
T = 1  # 周期

# 计算不同项数的傅里叶级数逼近
y1 = fourier_series_square(t, T, 1)
y3 = fourier_series_square(t, T, 3)
y5 = fourier_series_square(t, T, 5)
y50 = fourier_series_square(t, T, 50)
y_actual = square_wave(t, T)

# 绘制结果
plt.figure(figsize=(12, 10))

plt.subplot(3, 1, 1)
plt.plot(t, y_actual, 'k-', label='实际方波')
plt.plot(t, y1, 'r--', label='1项')
plt.legend()
plt.title('方波信号的傅里叶级数逼近 (1项)')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(t, y_actual, 'k-', label='实际方波')
plt.plot(t, y3, 'g--', label='3项')
plt.legend()
plt.title('方波信号的傅里叶级数逼近 (3项)')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(t, y_actual, 'k-', label='实际方波')
plt.plot(t, y50, 'b--', label='50项')
plt.legend()
plt.title('方波信号的傅里叶级数逼近 (50项)')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 生成时间序列
t = np.linspace(-np.pi, np.pi, 1000)

# 生成奇函数(正弦波)和偶函数(余弦波)
y_sin = np.sin(t)  # 奇函数
y_cos = np.cos(t)  # 偶函数

# 计算正弦波的傅里叶级数(只应包含正弦项)
def fourier_series_sin(t, n_terms):
    result = 0
    for n in range(1, n_terms+1):
        # 对于sin(t)，只有n=1时有系数，其他n的系数为0
        if n == 1:
            result += 2/np.pi * np.sin(n*t)
    return result

# 计算余弦波的傅里叶级数(只应包含余弦项和直流分量)
def fourier_series_cos(t, n_terms):
    result = 2/np.pi  # 直流分量
    for n in range(1, n_terms+1):
        # 特殊处理n=1的情况，避免除以零
        if n == 1:
            # cos(t)的傅里叶级数中n=1的系数为 t*sin(t)/2，需要单独处理
            result += 0.5 * t * np.sin(t)
        else:
            result += 4/(np.pi*(1 - n**2)) * np.cos(n*t)
    return result

# 计算逼近结果
y_sin_approx = fourier_series_sin(t, 5)
y_cos_approx = fourier_series_cos(t, 5)

# 绘制结果
plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)
plt.plot(t, y_sin, 'k-', label='sin(t)')
plt.plot(t, y_sin_approx, 'r--', label='傅里叶级数逼近')
plt.legend()
plt.title('奇函数(sin(t))及其傅里叶级数逼近')
plt.xlabel('时间 (rad)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 2, 2)
plt.plot(t, y_cos, 'k-', label='cos(t)')
plt.plot(t, y_cos_approx, 'b--', label='傅里叶级数逼近')
plt.legend()
plt.title('偶函数(cos(t))及其傅里叶级数逼近')
plt.xlabel('时间 (rad)')
plt.ylabel('振幅')
plt.grid(True)

# 验证奇偶性
plt.subplot(2, 2, 3)
plt.plot(t, y_sin, 'k-', label='sin(t)')
plt.plot(-t, -y_sin, 'r--', label='-sin(-t)')
plt.legend()
plt.title('奇函数验证: sin(-t) = -sin(t)')
plt.xlabel('时间 (rad)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 2, 4)
plt.plot(t, y_cos, 'k-', label='cos(t)')
plt.plot(-t, y_cos, 'b--', label='cos(-t)')
plt.legend()
plt.title('偶函数验证: cos(-t) = cos(t)')
plt.xlabel('时间 (rad)')
plt.ylabel('振幅')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 生成周期方波信号
def square_wave(t, T):
    return np.where(np.mod(t, T) < T / 2, 1, -1)

# 计算指数形式傅里叶系数
def fourier_coeff_exp(f, T, n_terms):
    """计算指数形式傅里叶系数"""
    t = np.linspace(-T / 2, T / 2, 1000)  # 一个周期
    f_t = f(t, T)
    c_n = np.zeros(n_terms * 2 + 1, dtype=complex)  # 存储系数

    for n in range(-n_terms, n_terms + 1):
        exponent = -1j * 2 * np.pi * n * t / T
        integrand = f_t * np.exp(exponent)
        # 使用np.trapezoid替代np.trapz
        c_n[n + n_terms] = (1 / T) * np.trapezoid(integrand, t)  # 数值积分

    return c_n

# 从指数系数重建信号
def reconstruct_from_exp_coeffs(t, T, c_n, n_terms):
    """使用指数形式傅里叶系数重建信号"""
    result = 0
    for n in range(-n_terms, n_terms + 1):
        exponent = 1j * 2 * np.pi * n * t / T
        result += c_n[n + n_terms] * np.exp(exponent)
    return np.real(result)  # 取实部

# 时间范围
t = np.linspace(-2, 2, 1000)
T = 1  # 周期
n_terms = 10  # 傅里叶系数项数

# 计算指数形式傅里叶系数
c_n = fourier_coeff_exp(square_wave, T, n_terms)

# 重建信号
y_reconstructed = reconstruct_from_exp_coeffs(t, T, c_n, n_terms)
y_actual = square_wave(t, T)

# 绘制结果
plt.figure(figsize=(12, 8))

plt.subplot(2, 1, 1)
plt.plot(t, y_actual, 'k-', label='实际方波')
plt.plot(t, y_reconstructed, 'r--', label='指数形式傅里叶级数重建')
plt.legend()
plt.title('使用指数形式傅里叶级数重建方波信号')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

# 绘制傅里叶系数的幅度谱
n_values = np.arange(-n_terms, n_terms + 1)
# 修正plt.stem的参数，只传递必要的参数
plt.subplot(2, 1, 2)
plt.stem(n_values, np.abs(c_n), linefmt='r-', markerfmt='ro', basefmt='r-')
plt.title('指数形式傅里叶系数的幅度谱')
plt.xlabel('谐波次数 n')
plt.ylabel('|c_n|')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 生成周期矩形脉冲信号
def periodic_rect_pulse(t, T, tau):
    """
    生成周期矩形脉冲信号
    t: 时间数组
    T: 周期
    tau: 脉冲宽度
    """
    return np.where(np.mod(t, T) < tau / 2, 1, 0)

# 计算周期信号的傅里叶系数(三角形式)
def fourier_coeff_trig(f, T, tau, n_terms):
    """计算三角形式傅里叶系数"""
    t = np.linspace(-T / 2, T / 2, 1000)  # 一个周期
    f_t = f(t, T, tau)

    # 直流分量 a0 - 使用 np.trapezoid 替代 np.trapz
    a0 = (2 / T) * np.trapezoid(f_t, t)

    a_n = np.zeros(n_terms + 1)  # 余弦项系数
    b_n = np.zeros(n_terms + 1)  # 正弦项系数

    for n in range(1, n_terms + 1):
        w0 = 2 * np.pi / T
        exponent_cos = n * w0 * t
        exponent_sin = n * w0 * t

        integrand_cos = f_t * np.cos(exponent_cos)
        integrand_sin = f_t * np.sin(exponent_sin)

        # 使用 np.trapezoid 替代 np.trapz
        a_n[n] = (2 / T) * np.trapezoid(integrand_cos, t)
        b_n[n] = (2 / T) * np.trapezoid(integrand_sin, t)

    return a0, a_n, b_n

# 时间范围
t = np.linspace(-5, 5, 1000)
T = 2  # 周期
tau = 0.5  # 脉冲宽度
n_terms = 20  # 傅里叶系数项数

# 生成周期矩形脉冲
y_rect = periodic_rect_pulse(t, T, tau)

# 计算傅里叶系数
a0, a_n, b_n = fourier_coeff_trig(periodic_rect_pulse, T, tau, n_terms)

# 重建信号
w0 = 2 * np.pi / T
t_reconstruct = np.linspace(-5, 5, 1000)
y_reconstructed = a0 / 2  # 直流分量

for n in range(1, n_terms + 1):
    y_reconstructed += a_n[n] * np.cos(n * w0 * t_reconstruct) + b_n[n] * np.sin(n * w0 * t_reconstruct)

# 计算幅度谱
amplitude_spectrum = np.sqrt(a_n ** 2 + b_n ** 2)
frequency = np.arange(0, n_terms + 1) * w0 / (2 * np.pi)  # 转换为Hz

# 绘制结果
plt.figure(figsize=(12, 10))

plt.subplot(2, 1, 1)
plt.plot(t, y_rect, 'k-', label=f'周期矩形脉冲 (T={T}, τ={tau})')
plt.plot(t_reconstruct, y_reconstructed, 'r--', label='傅里叶级数重建')
plt.legend()
plt.title('周期矩形脉冲信号及其傅里叶级数重建')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 1, 2)
# 修正 plt.stem 的参数，使用关键字参数
plt.stem(frequency, amplitude_spectrum, linefmt='b-', markerfmt='bo', basefmt='b-')
plt.title('周期矩形脉冲的幅度谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 定义单位冲激函数(近似)
def delta_function(t, epsilon=1e-3):
    """近似单位冲激函数"""
    return np.where(np.abs(t) < epsilon/2, 1/epsilon, 0)

# 定义矩形脉冲函数
def rect_pulse(t, tau):
    """矩形脉冲函数"""
    return np.where(np.abs(t) < tau/2, 1, 0)

# 计算傅里叶变换(数值方法)
def fourier_transform(f, t):
    """数值计算傅里叶变换"""
    N = len(t)
    dt = t[1] - t[0]
    w = np.fft.fftfreq(N, d=dt) * 2 * np.pi  # 角频率
    F = np.fft.fft(f) * dt  # 傅里叶变换
    return w, F

# 时间范围
t = np.linspace(-10, 10, 1000)
dt = t[1] - t[0]

# 计算单位冲激函数的傅里叶变换
delta = delta_function(t)
w_delta, F_delta = fourier_transform(delta, t)

# 计算矩形脉冲的傅里叶变换
tau = 2
rect = rect_pulse(t, tau)
w_rect, F_rect = fourier_transform(rect, t)

# 绘制结果
plt.figure(figsize=(12, 10))

plt.subplot(2, 2, 1)
plt.plot(t, delta, 'k-')
plt.title('单位冲激函数(近似)')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 2, 2)
plt.plot(w_delta, np.abs(F_delta), 'r-')
plt.title('单位冲激函数的傅里叶变换(幅度谱)')
plt.xlabel('角频率 (rad/s)')
plt.ylabel('|F(jω)|')
plt.grid(True)

plt.subplot(2, 2, 3)
plt.plot(t, rect, 'k-')
plt.title(f'矩形脉冲 (τ={tau}s)')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 2, 4)
plt.plot(w_rect, np.abs(F_rect), 'b-')
plt.title('矩形脉冲的傅里叶变换(幅度谱)')
plt.xlabel('角频率 (rad/s)')
plt.ylabel('|F(jω)|')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 定义矩形脉冲函数(偶函数)
def rect_pulse(t, tau):
    """矩形脉冲函数"""
    return np.where(np.abs(t) < tau/2, 1, 0)

# 定义符号函数(奇函数)
def sign_function(t):
    """符号函数"""
    return np.where(t > 0, 1, np.where(t < 0, -1, 0))

# 计算傅里叶变换
def fourier_transform(f, t):
    """数值计算傅里叶变换"""
    N = len(t)
    dt = t[1] - t[0]
    w = np.fft.fftfreq(N, d=dt) * 2 * np.pi  # 角频率
    F = np.fft.fft(f) * dt  # 傅里叶变换
    return w, F

# 时间范围
t = np.linspace(-10, 10, 1000)
dt = t[1] - t[0]

# 矩形脉冲(偶函数)
tau = 2
rect = rect_pulse(t, tau)
w_rect, F_rect = fourier_transform(rect, t)

# 符号函数(奇函数)
sign = sign_function(t)
w_sign, F_sign = fourier_transform(sign, t)

# 绘制结果
plt.figure(figsize=(12, 10))

# 矩形脉冲及其傅里叶变换
plt.subplot(2, 2, 1)
plt.plot(t, rect, 'k-')
plt.title('矩形脉冲(偶函数)')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 2, 2)
plt.plot(w_rect, np.real(F_rect), 'b-', label='实部')
plt.plot(w_rect, np.imag(F_rect), 'r--', label='虚部')
plt.title('矩形脉冲的傅里叶变换')
plt.xlabel('角频率 (rad/s)')
plt.ylabel('F(jω)')
plt.legend()
plt.grid(True)

# 符号函数及其傅里叶变换
plt.subplot(2, 2, 3)
plt.plot(t, sign, 'k-')
plt.title('符号函数(奇函数)')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 2, 4)
plt.plot(w_sign, np.real(F_sign), 'b-', label='实部')
plt.plot(w_sign, np.imag(F_sign), 'r--', label='虚部')
plt.title('符号函数的傅里叶变换')
plt.xlabel('角频率 (rad/s)')
plt.ylabel('F(jω)')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 定义高斯函数
def gaussian(t, sigma):
    """高斯函数"""
    return np.exp(-(t**2)/(2*sigma**2))

# 计算傅里叶变换
def fourier_transform(f, t):
    """数值计算傅里叶变换"""
    N = len(t)
    dt = t[1] - t[0]
    w = np.fft.fftfreq(N, d=dt) * 2 * np.pi  # 角频率
    F = np.fft.fft(f) * dt  # 傅里叶变换
    return w, F

# 时间范围
t = np.linspace(-10, 10, 1000)
dt = t[1] - t[0]

# 高斯函数(σ=1)
sigma1 = 1
g1 = gaussian(t, sigma1)
w1, F1 = fourier_transform(g1, t)

# 尺度变换(σ=2，时域扩展)
sigma2 = 2
g2 = gaussian(t, sigma2)
w2, F2 = fourier_transform(g2, t)

# 时移特性(σ=1，t0=2)
t0 = 2
g_shift = gaussian(t-t0, sigma1)
w_shift, F_shift = fourier_transform(g_shift, t)

# 绘制结果
plt.figure(figsize=(12, 10))

# 尺度变换
plt.subplot(2, 2, 1)
plt.plot(t, g1, 'b-', label='σ=1')
plt.plot(t, g2, 'r-', label='σ=2')
plt.legend()
plt.title('高斯函数的尺度变换')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 2, 2)
plt.plot(w1, np.abs(F1), 'b-', label='σ=1')
plt.plot(w2, np.abs(F2), 'r-', label='σ=2')
plt.legend()
plt.title('高斯函数傅里叶变换的尺度特性')
plt.xlabel('角频率 (rad/s)')
plt.ylabel('|F(jω)|')
plt.grid(True)

# 时移特性
plt.subplot(2, 2, 3)
plt.plot(t, g1, 'b-', label='原信号')
plt.plot(t, g_shift, 'r-', label=f'时移t0={t0}s')
plt.legend()
plt.title('高斯函数的时移')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 2, 4)
plt.plot(w1, np.abs(F1), 'b-', label='原信号幅度谱')
plt.plot(w_shift, np.abs(F_shift), 'r-', label='时移信号幅度谱')
plt.legend()
plt.title('高斯函数时移后的幅度谱')
plt.xlabel('角频率 (rad/s)')
plt.ylabel('|F(jω)|')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 定义矩形脉冲函数
def rect_pulse(t, tau):
    """矩形脉冲函数"""
    return np.where(np.abs(t) < tau/2, 1, 0)

# 计算傅里叶变换
def fourier_transform(f, t):
    """数值计算傅里叶变换"""
    N = len(t)
    dt = t[1] - t[0]
    w = np.fft.fftfreq(N, d=dt) * 2 * np.pi  # 角频率
    F = np.fft.fft(f) * dt  # 傅里叶变换
    return w, F

# 时间范围
t = np.linspace(-10, 10, 1000)
dt = t[1] - t[0]

# 不同宽度的矩形脉冲
tau1 = 1
tau2 = 2
tau3 = 4

rect1 = rect_pulse(t, tau1)
rect2 = rect_pulse(t, tau2)
rect3 = rect_pulse(t, tau3)

# 计算傅里叶变换
w1, F1 = fourier_transform(rect1, t)
w2, F2 = fourier_transform(rect2, t)
w3, F3 = fourier_transform(rect3, t)

# 计算能量谱密度
ESD1 = np.abs(F1)**2
ESD2 = np.abs(F2)** 2
ESD3 = np.abs(F3)**2

# 计算信号总能量(时域和频域)
energy_time1 = np.sum(rect1**2) * dt
energy_freq1 = np.sum(ESD1) * (w1[1] - w1[0]) / (2*np.pi)

energy_time2 = np.sum(rect2**2) * dt
energy_freq2 = np.sum(ESD2) * (w2[1] - w2[0]) / (2*np.pi)

energy_time3 = np.sum(rect3**2) * dt
energy_freq3 = np.sum(ESD3) * (w3[1] - w3[0]) / (2*np.pi)

print(f"τ={tau1}s 时，时域能量 = {energy_time1:.6f}, 频域能量 = {energy_freq1:.6f}")
print(f"τ={tau2}s 时，时域能量 = {energy_time2:.6f}, 频域能量 = {energy_freq2:.6f}")
print(f"τ={tau3}s 时，时域能量 = {energy_time3:.6f}, 频域能量 = {energy_freq3:.6f}")

# 绘制能量谱
plt.figure(figsize=(10, 6))
plt.plot(w1, ESD1, 'b-', label=f'τ={tau1}s')
plt.plot(w2, ESD2, 'g-', label=f'τ={tau2}s')
plt.plot(w3, ESD3, 'r-', label=f'τ={tau3}s')
plt.title('矩形脉冲信号的能量谱密度(ESD)')
plt.xlabel('角频率 (rad/s)')
plt.ylabel('能量谱密度 |F(jω)|^2')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 定义正弦信号
def sine_signal(t, f0):
    """正弦信号"""
    return np.sin(2 * np.pi * f0 * t)

# 计算傅里叶变换(使用FFT近似)
def fourier_transform(f, t):
    """数值计算傅里叶变换"""
    N = len(t)
    dt = t[1] - t[0]
    f_fft = np.fft.fft(f)
    f_fft_shifted = np.fft.fftshift(f_fft)
    freq = np.fft.fftfreq(N, d=dt)
    freq_shifted = np.fft.fftshift(freq)
    return freq_shifted, f_fft_shifted * dt

# 时间范围
t = np.linspace(-5, 5, 10000)
dt = t[1] - t[0]

# 正弦信号(f0=1Hz)
f0 = 1
sine = sine_signal(t, f0)

# 计算傅里叶变换
freq, f_fft = fourier_transform(sine, t)

# 绘制结果
plt.figure(figsize=(12, 8))

plt.subplot(2, 1, 1)
plt.plot(t, sine, 'k-')
plt.title(f'正弦信号 sin(2π*{f0}t)')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(freq, np.abs(f_fft), 'r-')
plt.title('正弦信号的傅里叶变换(幅度谱)')
plt.xlabel('频率 (Hz)')
plt.ylabel('|F(jω)|')
plt.grid(True)

# 标记理论上的冲激位置
plt.axvline(x=-f0, color='g', linestyle='--', label=f'理论冲激位置: ±{f0}Hz')
plt.axvline(x=f0, color='g', linestyle='--')
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 定义阶跃信号
def step_signal(t):
    """单位阶跃信号"""
    return np.where(t >= 0, 1, 0)

# 设计理想低通滤波器
def ideal_lowpass_filter(w, w_c):
    """理想低通滤波器频率响应"""
    return np.where(np.abs(w) < w_c, 1, 0)

# 计算傅里叶变换
def fourier_transform(f, t):
    """数值计算傅里叶变换"""
    N = len(t)
    dt = t[1] - t[0]
    w = np.fft.fftfreq(N, d=dt) * 2 * np.pi  # 角频率
    F = np.fft.fft(f) * dt  # 傅里叶变换
    return w, F

# 计算逆傅里叶变换
def inverse_fourier_transform(F, w):
    """数值计算逆傅里叶变换"""
    N = len(F)
    dw = w[1] - w[0]
    f = np.fft.ifft(F) * N * dw / (2 * np.pi)  # 逆傅里叶变换
    return f

# 时间范围
t = np.linspace(-10, 10, 10000)
dt = t[1] - t[0]

# 阶跃信号
step = step_signal(t)

# 计算阶跃信号的傅里叶变换
w_step, F_step = fourier_transform(step, t)

# 设计不同截止频率的理想低通滤波器
w_c1 = 1  # 低截止频率
w_c2 = 10  # 高截止频率

H1 = ideal_lowpass_filter(w_step, w_c1)
H2 = ideal_lowpass_filter(w_step, w_c2)

# 计算滤波后的信号
F_filtered1 = F_step * H1
F_filtered2 = F_step * H2

# 逆傅里叶变换得到时域信号
y_filtered1 = inverse_fourier_transform(F_filtered1, w_step)
y_filtered2 = inverse_fourier_transform(F_filtered2, w_step)

# 绘制结果
plt.figure(figsize=(12, 10))

plt.subplot(3, 1, 1)
plt.plot(t, step, 'k-')
plt.title('单位阶跃信号')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(w_step, np.abs(F_step), 'b-', label='阶跃信号频谱')
plt.plot(w_step, H1, 'r--', label=f'低通滤波器 (ωc={w_c1}rad/s)')
plt.plot(w_step, H2, 'g--', label=f'低通滤波器 (ωc={w_c2}rad/s)')
plt.title('阶跃信号频谱和滤波器频率响应')
plt.xlabel('角频率 (rad/s)')
plt.ylabel('幅度')
plt.legend()
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(t, step, 'k--', label='原始阶跃信号')
plt.plot(t, np.real(y_filtered1), 'r-', label=f'滤波后(ωc={w_c1})')
plt.plot(t, np.real(y_filtered2), 'b-', label=f'滤波后(ωc={w_c2})')
plt.title('理想低通滤波器对阶跃信号的响应')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False


# 定义带限信号(正弦波之和)
def bandlimited_signal(t, f1, f2):
    """带限信号: sin(2πf1t) + sin(2πf2t)"""
    return np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)


# 信号取样
def sample_signal(f, f1, f2, t, fs):
    """对信号进行取样"""
    Ts = 1 / fs  # 取样周期
    samples = f(t, f1, f2)
    t_samples = t  # 原始时间数组

    # 计算取样索引
    sample_indices = np.round(t / Ts).astype(int)
    unique_indices, inverse = np.unique(sample_indices, return_inverse=True)
    valid_indices = unique_indices[unique_indices < len(t)]

    return t[valid_indices], samples[valid_indices]


# 从取样信号重建信号(零阶保持)
def reconstruct_zero_order(t, t_samples, samples, fs):
    """使用零阶保持法从取样信号重建"""
    Ts = 1 / fs
    reconstructed = np.zeros_like(t)
    for i, ts in enumerate(t_samples):
        idx = np.where((t >= ts) & (t < ts + Ts))
        if idx[0].size > 0:  # 确保有索引可赋值
            reconstructed[idx] = samples[i]
    return reconstructed


# 时间范围
t = np.linspace(0, 2, 1000)
dt = t[1] - t[0]  # 时间间隔

# 带限信号参数
f1 = 1
f2 = 5  # 最高频率

# 不同的取样频率
fs_nyquist = 2 * f2  # 奈奎斯特频率(10Hz)
fs_undersample = f2  # 欠取样(5Hz)
fs_oversample = 10 * f2  # 过取样(50Hz)

# 取样
t_samples_nyquist, samples_nyquist = sample_signal(bandlimited_signal, f1, f2, t, fs_nyquist)
t_samples_undersample, samples_undersample = sample_signal(bandlimited_signal, f1, f2, t, fs_undersample)
t_samples_oversample, samples_oversample = sample_signal(bandlimited_signal, f1, f2, t, fs_oversample)

# 重建信号
reconstructed_nyquist = reconstruct_zero_order(t, t_samples_nyquist, samples_nyquist, fs_nyquist)
reconstructed_undersample = reconstruct_zero_order(t, t_samples_undersample, samples_undersample, fs_undersample)
reconstructed_oversample = reconstruct_zero_order(t, t_samples_oversample, samples_oversample, fs_oversample)

# 绘制结果
plt.figure(figsize=(12, 10))

plt.subplot(3, 1, 1)
plt.plot(t, bandlimited_signal(t, f1, f2), 'k-', label='原始带限信号')
# 修正plt.stem的参数，使用关键字参数
plt.stem(t_samples_nyquist, samples_nyquist, linefmt='r-', markerfmt='ro', basefmt='r-',
         label=f'取样点 (fs={fs_nyquist}Hz)')
plt.plot(t, reconstructed_nyquist, 'b--', label='奈奎斯特频率重建')
plt.legend()
plt.title('奈奎斯特频率取样和重建')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(t, bandlimited_signal(t, f1, f2), 'k-', label='原始带限信号')
plt.stem(t_samples_undersample, samples_undersample, linefmt='r-', markerfmt='ro', basefmt='r-',
         label=f'取样点 (fs={fs_undersample}Hz)')
plt.plot(t, reconstructed_undersample, 'g--', label='欠取样重建')
plt.legend()
plt.title('欠取样(fs < 2fm)时的信号失真')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(t, bandlimited_signal(t, f1, f2), 'k-', label='原始带限信号')
plt.stem(t_samples_oversample, samples_oversample, linefmt='r-', markerfmt='ro', basefmt='r-',
         label=f'取样点 (fs={fs_oversample}Hz)')
plt.plot(t, reconstructed_oversample, 'm--', label='过取样重建')
plt.legend()
plt.title('过取样(fs > 2fm)时的信号重建')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 定义周期矩形序列
def periodic_rect_sequence(n, N, M):
    """
    周期矩形序列
    n: 时间索引
    N: 周期
    M: 矩形脉冲宽度(M < N)
    """
    return np.where(np.mod(n, N) < M, 1, 0)

# 计算离散傅里叶级数(DFS)
def dfs(x, N):
    """计算离散傅里叶级数"""
    n = np.arange(N)
    k = np.arange(N)
    X = np.zeros(N, dtype=complex)

    for i in range(N):
        X[i] = np.sum(x * np.exp(-1j * 2 * np.pi * k[i] * n / N))

    return X

# 从DFS系数重建序列
def idfs(X, N):
    """计算逆离散傅里叶级数"""
    n = np.arange(N)
    k = np.arange(N)
    x = np.zeros(N, dtype=complex)

    for i in range(N):
        x[i] = (1 / N) * np.sum(X * np.exp(1j * 2 * np.pi * k[i] * n / N))

    return x

# 序列长度和周期
N = 16  # 周期
M = 4  # 矩形脉冲宽度
n = np.arange(3 * N)  # 绘制3个周期

# 生成周期矩形序列
x = periodic_rect_sequence(n, N, M)

# 计算DFS
X = dfs(x[:N], N)

# 从DFS重建序列
x_reconstructed = idfs(X, N)

# 绘制结果
plt.figure(figsize=(12, 8))

plt.subplot(2, 1, 1)
# 修改后的stem函数调用
plt.stem(n, x, linefmt='b-', markerfmt='bo', basefmt='b-')
plt.title(f'周期矩形序列 (N={N}, M={M})')
plt.xlabel('时间索引 n')
plt.ylabel('x(n)')
plt.grid(True)

plt.subplot(2, 1, 2)
# 修改后的stem函数调用
plt.stem(np.arange(N), np.abs(X), linefmt='r-', markerfmt='ro', basefmt='r-')
plt.title('离散傅里叶级数(DFS)幅度谱')
plt.xlabel('频率索引 k')
plt.ylabel('|X(k)|')
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 生成测试信号(正弦波+噪声)
def generate_test_signal(t, f1, f2, snr=20):
    """生成测试信号: 两个正弦波+噪声"""
    signal = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)

    # 计算信号功率
    signal_power = np.mean(signal ** 2)

    # 计算噪声功率以达到指定SNR
    noise_power = signal_power / (10 ** (snr / 10))

    # 生成噪声
    noise = np.random.normal(0, np.sqrt(noise_power), len(t))

    # 带噪信号
    noisy_signal = signal + noise

    return noisy_signal, signal, noise


# 使用FFT计算频谱
def calculate_spectrum(x, fs):
    """使用FFT计算信号频谱"""
    N = len(x)
    X = np.fft.fft(x)
    X_mag = np.abs(X) / N  # 幅度谱
    freq = np.fft.fftfreq(N, d=1 / fs)  # 频率轴
    return freq, X_mag


# 时间参数
fs = 1000  # 取样频率(Hz)
t = np.linspace(0, 1, fs)  # 1秒信号

# 信号频率
f1 = 50
f2 = 150

# 生成测试信号
noisy_signal, clean_signal, noise = generate_test_signal(t, f1, f2, snr=15)

# 计算频谱
freq_noisy, X_mag_noisy = calculate_spectrum(noisy_signal, fs)
freq_clean, X_mag_clean = calculate_spectrum(clean_signal, fs)

# 只考虑正频率部分
pos_freq_noisy = freq_noisy[freq_noisy >= 0]
X_mag_noisy_pos = X_mag_noisy[freq_noisy >= 0]

pos_freq_clean = freq_clean[freq_clean >= 0]
X_mag_clean_pos = X_mag_clean[freq_clean >= 0]

# 绘制结果
plt.figure(figsize=(12, 10))

plt.subplot(2, 1, 1)
plt.plot(t, clean_signal, 'b-', label='干净信号')
plt.plot(t, noisy_signal, 'r--', label='带噪信号')
plt.legend()
plt.title('测试信号(干净信号和带噪信号)')
plt.xlabel('时间 (s)')
plt.ylabel('振幅')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(pos_freq_clean, X_mag_clean_pos, 'b-', label='干净信号频谱')
plt.plot(pos_freq_noisy, X_mag_noisy_pos, 'r--', label='带噪信号频谱')
plt.axvline(x=f1, color='g', linestyle='--', label=f'信号频率: {f1}Hz')
plt.axvline(x=f2, color='g', linestyle='--', label=f'信号频率: {f2}Hz')
plt.legend()
plt.title('使用FFT计算的信号频谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.grid(True)

plt.tight_layout()
plt.show()