import numpy as np from math import * def Pnm(Phi, Degree): P = np.zeros([Degree + 2, Degree + 2]) # 跨阶次正规化勒让德系数 P[1][1] = 1 P[2][1] = sin(Phi) * 3 ** 0.5 P[2][2] = sqrt(3 * (1 - sin(Phi) ** 2)) for j in range(1, 3): for i in range(3, Degree + 2): l = i - 1 m = j - 1 a = sqrt((4 * l ** 2 - 1) / (l ** 2 - m ** 2)) b = sqrt((2 * l + 1) / (2 * l - 3)) * sqrt(((l - 1) ** 2 - m ** 2) / (l ** 2 - m ** 2)) P[i][j] = a * sin(Phi) * P[i - 1][j] - b * P[i - 2][j] for j in range(3, Degree + 1): for i in range(j, j + 2): l = i - 1 m = j - 1 if (m == 2): beta = sqrt(2 * (2 * l + 1) * (l + m - 2) * (l + m - 3) / (2 * l - 3) / (l + m) / (l + m - 1)) gama = sqrt(2 * (l - m + 1) * (l - m + 2) / (l + m) / (l + m - 1)) else: beta = sqrt((2 * l + 1) * (l + m - 2) * (l + m - 3) / (2 * l - 3) / (l + m) / (l + m - 1)) gama = sqrt((l - m + 1) * (l - m + 2) / (l + m) / (l + m - 1)) P[i][j] = beta * P[i - 2][j - 2] - gama * P[i][j - 2] if ((j + 2) < Degree + 2): for i in range(j + 2, Degree + 2): l = i - 1 m = j - 1 alpha = sqrt((2 * l + 1) * (l - m) * (l - m - 1) / (2 * l - 3) / (l + m) / (l + m - 1)) if (m == 2): beta = sqrt(2 * (2 * l + 1) * (l + m - 2) * (l + m - 3) / (2 * l - 3) / (l + m) / (l + m - 1)) gama = sqrt(2 * (l - m + 1) * (l - m + 2) / (l + m) / (l + m - 1)) else: beta = sqrt((2 * l + 1) * (l + m - 2) * (l + m - 3) / (2 * l - 3) / (l + m) / (l + m - 1)) gama = sqrt((l - m + 1) * (l - m + 2) / (l + m) / (l + m - 1)) P[i][j] = alpha * P[i - 2][j] + beta * P[i - 2][j - 2] - gama * P[i][j - 2] l = Degree m = Degree beta = sqrt((2 * l + 1) * (l + m - 2) * (l + m - 3) / (2 * l - 3) / (l + m) / (l + m - 1)) gama = sqrt((l - m + 1) * (l - m + 2) / (l + m) / (l + m - 1)) P[l + 1][m + 1] = beta * P[l + 1 - 2][m + 1 - 2] - gama * P[l + 1][m + 1 - 2] return P def P_final(theta, n, m, Degree=360): Phi = pi / 2 - theta res = Pnm(Phi, Degree) return res a = P_final(radians(58), 360, 360) print(a)

for n in range(Degree): for m in range(n, Degree): if n == m: P[n, m] = sqrt((2 - 1) / 2 * factorial(n) / (4 * pi * factorial(n))) * cos(m * Phi) else: P[n, m] = sqrt((2 * n + 1) / (2 * n * (n + 1)) *

相关问题

优化这段import numpy as np import matplotlib.pyplot as plt %config InlineBackend.figure_format='retina' def generate_signal(t_vec, A, phi, noise, freq): Omega = 2*np.pi*freq return A * np.sin(Omega*t_vec + phi) + noise * (2*np.random.random def lock_in_measurement(signal, t_vec, ref_freq): Omega = 2*np.pi*ref_freq ref_0 = 2*np.sin(Omega*t_vec) ref_1 = 2*np.cos(Omega*t_vec) # signal_0 = signal * ref_0 signal_1 = signal * ref_1 # X = np.mean(signal_0) Y = np.mean(signal_1) # A = np.sqrt(X**2+Y**2) phi = np.arctan2(Y,X) print("A=", A, "phi=", phi) # t_vec = np.linspace(0, 0.2, 1001) A = 1 phi = np.pi noise = 0.2 ref_freq = 17.77777 # signal = generate_signal(t_vec, A, phi, noise, ref_freq) # lock_in_measurement(signal, t_vec, ref_freq)

import numpy as np
import matplotlib.pyplot as plt

%config InlineBackend.figure_format='retina'

def generate_signal(t_vec, A, phi, noise, freq):
Omega = 2*np.pi*freq
return A * np.sin(Omega*t_vec + phi) + noise * (2*np.random.random)

def lock_in_measurement(signal, t_vec, ref_freq):
Omega = 2*np.pi*ref_freq
ref_0 = 2*np.sin(Omega*t_vec)
ref_1 = 2*np.cos(Omega*t_vec)

signal_0 = signal * ref_0
signal_1 = signal * ref_1

X = np.mean(signal_0)
Y = np.mean(signal_1)

A = np.sqrt(X**2+Y**2)
phi = np.arctan2(Y,X)

print("A=", A, "phi=", phi)

t_vec = np.linspace(0, 0.2, 1001)
A = 1
phi = np.pi
noise = 0.2
ref_freq = 17.77777

signal = generate_signal(t_vec, A, phi, noise, ref_freq)
lock_in_measurement(signal, t_vec, ref_freq)

优化这段pythonimport numpy as np import matplotlib.pyplot as plt import math # 待测信号 freq = 17.77777 # 信号频率 t = np.linspace(0, 0.2, 1001) Omega =2 * np.pi * freq phi = np.pi A=1 x = A * np.sin(Omega * t + phi) # 加入噪声 noise = 0.2 * np.random.randn(len(t)) x_noi

se = x + noise # 绘制原始信号和加噪声后的信号 plt.figure(figsize=(10, 4)) plt.plot(t, x, label='Original Signal') plt.plot(t, x_noise, label='Signal with Noise') plt.legend() plt.xlabel('Time (s)') plt.ylabel('Amplitude') plt.show() # 进行傅里叶变换 fft_x_noisese = np.fft.fft(x_noise) freqs = np.fft.fftfreq(len(x_noise)) # 绘制频谱图 plt.figure(figsize=(10, 4)) plt.plot(freqs, np.abs(fft_x_noisese)) plt.xlabel('Frequency (Hz)') plt.ylabel('Amplitude') plt.show()

优化建议：
1. 可以将一些常量提取出来，例如频率、噪声幅度等，避免在循环中重复计算。
2. 可以使用subplot函数将多张图放在同一张画布中展示，提高可视化效率。
3. 可以对频谱图进行对数变换，使其更容易观察信号的频域特征。
4. 可以对傅里叶变换结果进行归一化处理，使得频谱图的纵轴单位更易理解。