A:
I have a theta_a_lst with size 1×721, It also represents one cycle of 50Hz sin wave. I also have voltage waveform va_lst. Perform FFT on va_lst.
Q:
Got it. Since you have one cycle of a 50 Hz sine (theta_a_lst size 1×721) and a corresponding voltage waveform va_lst, here are concise MATLAB and Python options to compute the FFT, report amplitudes/phases, and plot up to Nyquist.
Assumptions
– theta_a_lst spans exactly one electrical cycle (0 to 2π or 0 to 360°) with 721 samples.
– va_lst is the same length as theta_a_lst.
– Effective sampling frequency Fs is derived from 50 Hz fundamental: Fs = 50 × N (samples per cycle).
– Frequency resolution df = Fs/N = 50 Hz (harmonic bins at 0, 50, 100, …).
这里的有效采样频率（effective sampling frequency）Fs很重要，theta_lst为一个50Hz基波周期，其中点数N越大，Fs越大，但是频率分辨率df却是不变的，因为Fs = 50xN, 而df = Fs/N, 所以df = 50Hz。
如果只给一个50Hz的基波所对应的序列，能分析出来的谐波肯定都是50Hz的倍数。
MATLAB version
function out = fft_va_from_theta(theta_a_lst, va_lst, varargin)
% fft_va_from_theta FFT of va given one-cycle theta vector (50 Hz fundamental)
% Usage:
% out = fft_va_from_theta(theta_a_lst, va_lst, 'FmaxPlot', 2500, 'DoPlot', true)
%
% Inputs:
% theta_a_lst : 1xN or Nx1 angles covering one cycle (rad or deg)
% va_lst : 1xN or Nx1 voltage samples corresponding to theta_a_lst
%
% Name-Value:
% 'F0' : fundamental frequency (default 50 Hz)
% 'Window' : window name or vector (default 'hann')
% 'DoPlot' : true/false (default true)
% 'FmaxPlot' : max frequency for plots (default min(2500, Fs/2))
%
% Output struct out:
% .Fs, .df, .N, .f1 (kmaxx1), .X1 (kmaxx1 complex)
% .A1 (single-sided amplitude), .Sxx1 (single-sided PSD, arbitrary units if no Ts)
% .DC, .peak_amp, .peak_freq
% .idx_harm (indices at integer harmonics), .A1_harm, .phase_harm
p = inputParser;
p.addParameter('F0', 50, @(x)isnumeric(x)&&isscalar(x)&&x>0);
p.addParameter('Window','hann');
p.addParameter('DoPlot', true, @(x)islogical(x)&&isscalar(x));
p.addParameter('FmaxPlot', 2500, @(x)isnumeric(x)&&isscalar(x)&&x>0);
p.parse(varargin{:});
F0 = p.Results.F0;
wopt = p.Results.Window;
DoPlot = p.Results.DoPlot;
% Ensure column vectors
theta = theta_a_lst(:);
v = va_lst(:);
N = numel(v);
assert(numel(theta)==N, 'theta_a_lst and va_lst must have same length');
% Sampling frequency from one-cycle assumption
Fs = F0 * N;
df = Fs / N; % equals F0
% Window
if ischar(wopt) || (isstring(wopt) && isscalar(wopt))
switch lower(string(wopt))
case "hann"
w = hann(N, 'periodic');
case "hamming"
w = hamming(N, 'periodic');
case "rect"
w = ones(N,1);
otherwise
error('Unsupported Window option.');
end
else
w = wopt(:);
assert(numel(w)==N, 'Custom window must match data length');
end
W0 = sum(w)/N; % coherent gain
Pw = mean(w.^2); % window power (for PSD if needed)
vw = v .* w;
% FFT
X = fft(vw);
Nfft = N;
f = (0:Nfft-1).' * (Fs/Nfft);
% Single-sided
if mod(Nfft,2)==0
kmax = Nfft/2 + 1;
else
kmax = (Nfft+1)/2;
end
f1 = f(1:kmax);
X1 = X(1:kmax);
% Amplitude spectrum (single-sided), window-corrected
A2 = abs(X1) / (N * W0);
A1 = A2;
if kmax > 2
A1(2:end-1) = 2*A1(2:end-1);
end
% PSD (single-sided). Units V^2/Hz if Ts known; here Ts=1/Fs.
Sxx2 = (abs(X1).^2) / (Fs * N * Pw);
Sxx1 = Sxx2;
if kmax > 2
Sxx1(2:end-1) = 2*Sxx1(2:end-1);
end
% Stats
DC = A1(1);
[peak_amp, idxp] = max(A1(2:end)); idxp = idxp + 1;
peak_freq = f1(idxp);
% Harmonic bins (0, F0, 2F0, ...)
harmonics = (0:floor(f1(end)/F0)).' * F0;
[~, idx_harm] = ismember(harmonics, f1); % exact because df=F0
A1_harm = A1(idx_harm);
phase_harm = angle(X1(idx_harm));
% Pack outputs
out = struct('Fs',Fs,'df',df,'N',N,'f1',f1,'X1',X1, ...
'A1',A1,'Sxx1',Sxx1,'DC',DC, ...
'peak_amp',peak_amp,'peak_freq',peak_freq, ...
'idx_harm',idx_harm,'A1_harm',A1_harm,'phase_harm',phase_harm, ...
'W0',W0,'Pw',Pw);
% Plots
if DoPlot
fmax_plot = min(2500, Fs/2);
idx_plot = f1 <= fmax_plot;
fig = figure('Name','va FFT (one-cycle, 50 Hz)','Color','w');
% scale 50%
pos = get(fig,'Position'); pos(3:4) = 0.5*pos(3:4); set(fig,'Position',pos);
tiledlayout(3,1);
nexttile;
plot((0:N-1)/Fs, v, 'k-'); grid on;
xlabel('Time (s)'); ylabel('v_a');
title('Time-domain (one cycle)');
nexttile;
stem(f1(idx_plot), A1(idx_plot), 'filled'); grid on;
xlim([0, fmax_plot]); xlabel('Frequency (Hz)'); ylabel('Amplitude (V)');
title(sprintf('Single-Sided Amplitude Spectrum (df = %.1f Hz)', df));
nexttile;
plot(f1(idx_plot), 10*log10(Sxx1(idx_plot)+eps), 'b-'); grid on;
xlim([0, fmax_plot]); xlabel('Frequency (Hz)'); ylabel('PSD (dBV^2/Hz)');
title('Single-Sided PSD');
end
end
Python version (NumPy)
import numpy as np
def fft_va_from_theta(theta_a_lst, va_lst, f0=50.0, window='hann'):
"""
FFT of va given one-cycle theta vector (50 Hz fundamental).
Returns dict with Fs, df, f1, X1, A1, phases, etc.
"""
theta = np.asarray(theta_a_lst).ravel()
v = np.asarray(va_lst).ravel()
assert theta.size == v.size
N = v.size
Fs = f0 * N
df = Fs / N # equals f0
if isinstance(window, str):
if window.lower() == 'hann':
w = np.hanning(N)
elif window.lower() == 'hamming':
w = np.hamming(N)
elif window.lower() == 'rect':
w = np.ones(N)
else:
raise ValueError('Unsupported window')
else:
w = np.asarray(window).ravel()
assert w.size == N
W0 = w.sum()/N
Pw = np.mean(w**2)
vw = v * w
X = np.fft.fft(vw)
Nfft = N
f = np.arange(Nfft) * (Fs/Nfft)
kmax = Nfft//2 + 1 if Nfft % 2 == 0 else (Nfft+1)//2
f1 = f[:kmax]
X1 = X[:kmax]
A2 = np.abs(X1) / (N * W0)
A1 = A2.copy()
if kmax > 2:
A1[1:-1] *= 2.0
Sxx2 = (np.abs(X1)**2) / (Fs * N * Pw)
Sxx1 = Sxx2.copy()
if kmax > 2:
Sxx1[1:-1] *= 2.0
DC = A1[0]
idxp = 1 + np.argmax(A1[1:])
peak_amp = A1[idxp]
peak_freq = f1[idxp]
harmonics = np.arange(int(f1[-1]//f0)+1) * f0
idx_harm = (harmonics / df).astype(int)
A1_harm = A1[idx_harm]
phase_harm = np.angle(X1[idx_harm])
return dict(Fs=Fs, df=df, N=N, f1=f1, X1=X1, A1=A1, Sxx1=Sxx1, DC=DC,
peak_amp=peak_amp, peak_freq=peak_freq, idx_harm=idx_harm,
A1_harm=A1_harm, phase_harm=phase_harm, W0=W0, Pw=Pw)
If you share the units/format of theta_a_lst (deg vs rad) or want a specific window, I can tailor the code further.