通过谐波数生成直方图
Generating a Histogram by Harmonic Number
我正在尝试在 GNU Octave 中创建一个程序来绘制一个直方图,显示修改后的正弦波的基波和谐波(SCR 调光器的输出,它由一个正弦波组成,该正弦波在部分通过波)。
我已经能够生成波形并执行 FFT 以获得一组频率与幅度点,但是我不确定如何将这些数据转换为适合生成直方图的 bin。
示例代码和下面我所追求的图像 - 感谢您的帮助!
clear();
vrms = 120;
freq = 60;
nCycles = 2;
level = 25;
vpeak = sqrt(2) * vrms;
sampleinterval = 0.00001;
num_harmonics = 10
disp("Start");
% Draw the waveform
x = 0 : sampleinterval : nCycles * 1 / freq; % time in sampleinterval increments
dimmed_wave = [];
undimmed_wave = [];
for i = 1 : columns(x)
rad_value = x(i) * 2 * pi * freq;
off_time = mod(rad_value, pi);
on_time = pi*(100-level)/100;
if (off_time < on_time)
dimmed_wave = [dimmed_wave, 0]; % in the dimmed period, value is zero
else
dimmed_wave = [dimmed_wave, sin(rad_value)]; % when not dimmed, value = sine
endif
undimmed_wave = [undimmed_wave, sin(rad_value)];
endfor
y = dimmed_wave * vpeak; % calculate instantaneous voltage
undimmed = undimmed_wave * vpeak;
subplot(2,1,1)
plot(x*1000, y, '-', x*1000, undimmed, '--');
xlabel ("Time (ms)");
ylabel ("Voltage");
% Fourier Transform to determine harmonics
subplot(2,1,2)
N = length(dimmed_wave); % number of points
fft_vals = abs(fftshift(fft(dimmed_wave))); % perform fft
frequency = [ -(ceil((N-1)/2):-1:1) ,0 ,(1:floor((N-1)/2)) ] * 1 / (N *sampleinterval);
plot(frequency, fft_vals);
axis([0,400]);
xlabel ("Frequency");
ylabel ("Amplitude");
你知道你的基频(基音),我们称它为F
。 2*F
是二次谐波,3*F
是三次谐波,等等。您想要在这些之间的中间设置直方图 bin 边缘:1.5*F
、2.5*F
等。
您的输入信号有两个周期,因此您的(整数)基频是 k=2
([=17= 处的值,图中的第一个峰值)。二次谐波在 k=4
,三次谐波在 k=6
,等等
因此您可以将 bin 边缘设置为 k = 1:2:end
。
一般来说,这会是 k = nCycles/2:nCycles:end
。
您可以根据我们计算的 bin 边缘计算您的条形图,如下所示:
fft_vals = abs(fft(dimmed_wave));
nHarmonics = 9;
edges = nCycles/2 + (0:nHarmonics)*nCycles;
H = cumsum(fft_vals);
H = diff(H(edges));
bar(1:nHarmonics,H);
我正在尝试在 GNU Octave 中创建一个程序来绘制一个直方图,显示修改后的正弦波的基波和谐波(SCR 调光器的输出,它由一个正弦波组成,该正弦波在部分通过波)。
我已经能够生成波形并执行 FFT 以获得一组频率与幅度点,但是我不确定如何将这些数据转换为适合生成直方图的 bin。
示例代码和下面我所追求的图像 - 感谢您的帮助!
clear();
vrms = 120;
freq = 60;
nCycles = 2;
level = 25;
vpeak = sqrt(2) * vrms;
sampleinterval = 0.00001;
num_harmonics = 10
disp("Start");
% Draw the waveform
x = 0 : sampleinterval : nCycles * 1 / freq; % time in sampleinterval increments
dimmed_wave = [];
undimmed_wave = [];
for i = 1 : columns(x)
rad_value = x(i) * 2 * pi * freq;
off_time = mod(rad_value, pi);
on_time = pi*(100-level)/100;
if (off_time < on_time)
dimmed_wave = [dimmed_wave, 0]; % in the dimmed period, value is zero
else
dimmed_wave = [dimmed_wave, sin(rad_value)]; % when not dimmed, value = sine
endif
undimmed_wave = [undimmed_wave, sin(rad_value)];
endfor
y = dimmed_wave * vpeak; % calculate instantaneous voltage
undimmed = undimmed_wave * vpeak;
subplot(2,1,1)
plot(x*1000, y, '-', x*1000, undimmed, '--');
xlabel ("Time (ms)");
ylabel ("Voltage");
% Fourier Transform to determine harmonics
subplot(2,1,2)
N = length(dimmed_wave); % number of points
fft_vals = abs(fftshift(fft(dimmed_wave))); % perform fft
frequency = [ -(ceil((N-1)/2):-1:1) ,0 ,(1:floor((N-1)/2)) ] * 1 / (N *sampleinterval);
plot(frequency, fft_vals);
axis([0,400]);
xlabel ("Frequency");
ylabel ("Amplitude");
你知道你的基频(基音),我们称它为F
。 2*F
是二次谐波,3*F
是三次谐波,等等。您想要在这些之间的中间设置直方图 bin 边缘:1.5*F
、2.5*F
等。
您的输入信号有两个周期,因此您的(整数)基频是 k=2
([=17= 处的值,图中的第一个峰值)。二次谐波在 k=4
,三次谐波在 k=6
,等等
因此您可以将 bin 边缘设置为 k = 1:2:end
。
一般来说,这会是 k = nCycles/2:nCycles:end
。
您可以根据我们计算的 bin 边缘计算您的条形图,如下所示:
fft_vals = abs(fft(dimmed_wave));
nHarmonics = 9;
edges = nCycles/2 + (0:nHarmonics)*nCycles;
H = cumsum(fft_vals);
H = diff(H(edges));
bar(1:nHarmonics,H);