Matlab 上的 DFT 代码无法按预期工作
DFT code on Matlab does not work as intended
我正在尝试在 Matlab 上实现基本的 DFT 算法。
我只是使用正弦波的相位和正交分量进行相位调制(增加频率 a.k.a 线性调频)。我确实将我的结果与 Matlab 的 fft 命令进行了比较。只要没有相位调制(纯正弦),我的代码就会给出相同的结果。每当我添加线性调频调制时,结果都会不同。例如,当我在载波周围使用具有一定带宽的线性调频时,预期结果应该是从载波频率开始的线性调频带宽的频率分布。但是,我也从载波频率开始向后获得该结果的副本。您可以在下面使用我的代码而无需修改任何内容。图 5 是我的结果,图 6 是预期结果。载波为 256 Hz,带有 10Hz 带宽的线性调频信号。你可以看到下面的代码。重要的部分是 for 循环,我在其中获取信号的 dft。你也可以在下面看到我的 dft 结果。
close all;
clear all;
%% signal generation
t = (0:0.0001:1); % 1 second window
f = 256; %freq of input signal in hertz
bw = 10; % bandwidth sweep of signal
phaseInput = 2*pi*t*bw.*t;
signalInput = sin(2*pi*f*t + phaseInput); %input signal
inphase = sin(2*pi*f*t).*cos(phaseInput); %inphase component
quadrature = cos(2*pi*f*t).*sin(phaseInput); %quadrature component
figure
plot(t,signalInput,'b',t,inphase,'g',t,quadrature,'r');
title('Input Signal');
xlabel('Time in seconds');
ylabel('Amplitude');
%% sampling signal previously generated
Fs = 1024; %sampling freq
Ts = (0:1/Fs:1);%sample times for 1 second window
sPhase = 2*pi*Ts*bw.*Ts;
sI = sin(2*pi*f*Ts).*cos(sPhase);
sQ = cos(2*pi*f*Ts).*sin(sPhase);
hold on;
plot(Ts,sI+sQ,'b*',Ts,sI,'g*',Ts,sQ,'r*');
fftSize = Fs; %Using all samples in dft
sampleIdx = (0:1:fftSize-1)';
sampledI = sI(1:fftSize)';
sampledQ = sQ(1:fftSize)';
figure;
plot(sampleIdx,sampledI,sampleIdx,sampledQ);
title('Sampled IQ Components');
%% DFT Calculation
dftI = zeros(fftSize,1);
dftQ = zeros(fftSize,1);
for w = 0:fftSize-1
%exp(-2*pi*w*t) = cos(2*pi*w*t) - i*sin(2*pi*w*t)
cI = cos(2*pi*w*sampleIdx/fftSize); %correlation cos
cQ = -sin(2*pi*w*sampleIdx/fftSize); %correlation sin
dftI(w+1) = sum(sampledI.*cI - sampledQ.*cQ); %
dftQ(w+1) = sum(sampledI.*cQ + sampledQ.*cI);
end;
figure;
plot(Fs*sampleIdx/fftSize,dftI);
title('DFT Inphase');
xlabel('Hertz');
figure
plot(Fs*sampleIdx/fftSize,dftQ);
title('DFT Quadrature');
xlabel('Hertz');
figure;
plot(Fs*sampleIdx/fftSize,sqrt(dftQ.^2+dftI.^2));
%% For comparison
sampledInput = sin(2*pi*f*Ts + sPhase);
Y = fft(sampledInput(1:1024),1024);
Pyy = Y.*conj(Y)/1024;
f = (0:1023);
figure;
plot(f,Pyy)
title('Power spectral density')
xlabel('Frequency (Hz)')
原因就在于两个不同的信号肯定会给你两个不同的频谱。查看下面的代码,你会发现你实际给出的dft算法的输入是sampledI+jsampledQ
。因此,您在这里所做的不仅仅是将原始信号分解为 In-phase and quadrature components
,而是在这里进行 Hilbert transform
——将 real
信号更改为 [=15] =]一个。
cI = cos(2*pi*w*sampleIdx/fftSize); %correlation cos
cQ = -sin(2*pi*w*sampleIdx/fftSize); %correlation sin
dftI(w+1) = sum(sampledI.*cI - sampledQ.*cQ); %
dftQ(w+1) = sum(sampledI.*cQ + sampledQ.*cI);
所以比较的 sampledInput
应该是 sampledInput = sI+1i*sQ;
.
我正在尝试在 Matlab 上实现基本的 DFT 算法。 我只是使用正弦波的相位和正交分量进行相位调制(增加频率 a.k.a 线性调频)。我确实将我的结果与 Matlab 的 fft 命令进行了比较。只要没有相位调制(纯正弦),我的代码就会给出相同的结果。每当我添加线性调频调制时,结果都会不同。例如,当我在载波周围使用具有一定带宽的线性调频时,预期结果应该是从载波频率开始的线性调频带宽的频率分布。但是,我也从载波频率开始向后获得该结果的副本。您可以在下面使用我的代码而无需修改任何内容。图 5 是我的结果,图 6 是预期结果。载波为 256 Hz,带有 10Hz 带宽的线性调频信号。你可以看到下面的代码。重要的部分是 for 循环,我在其中获取信号的 dft。你也可以在下面看到我的 dft 结果。
close all;
clear all;
%% signal generation
t = (0:0.0001:1); % 1 second window
f = 256; %freq of input signal in hertz
bw = 10; % bandwidth sweep of signal
phaseInput = 2*pi*t*bw.*t;
signalInput = sin(2*pi*f*t + phaseInput); %input signal
inphase = sin(2*pi*f*t).*cos(phaseInput); %inphase component
quadrature = cos(2*pi*f*t).*sin(phaseInput); %quadrature component
figure
plot(t,signalInput,'b',t,inphase,'g',t,quadrature,'r');
title('Input Signal');
xlabel('Time in seconds');
ylabel('Amplitude');
%% sampling signal previously generated
Fs = 1024; %sampling freq
Ts = (0:1/Fs:1);%sample times for 1 second window
sPhase = 2*pi*Ts*bw.*Ts;
sI = sin(2*pi*f*Ts).*cos(sPhase);
sQ = cos(2*pi*f*Ts).*sin(sPhase);
hold on;
plot(Ts,sI+sQ,'b*',Ts,sI,'g*',Ts,sQ,'r*');
fftSize = Fs; %Using all samples in dft
sampleIdx = (0:1:fftSize-1)';
sampledI = sI(1:fftSize)';
sampledQ = sQ(1:fftSize)';
figure;
plot(sampleIdx,sampledI,sampleIdx,sampledQ);
title('Sampled IQ Components');
%% DFT Calculation
dftI = zeros(fftSize,1);
dftQ = zeros(fftSize,1);
for w = 0:fftSize-1
%exp(-2*pi*w*t) = cos(2*pi*w*t) - i*sin(2*pi*w*t)
cI = cos(2*pi*w*sampleIdx/fftSize); %correlation cos
cQ = -sin(2*pi*w*sampleIdx/fftSize); %correlation sin
dftI(w+1) = sum(sampledI.*cI - sampledQ.*cQ); %
dftQ(w+1) = sum(sampledI.*cQ + sampledQ.*cI);
end;
figure;
plot(Fs*sampleIdx/fftSize,dftI);
title('DFT Inphase');
xlabel('Hertz');
figure
plot(Fs*sampleIdx/fftSize,dftQ);
title('DFT Quadrature');
xlabel('Hertz');
figure;
plot(Fs*sampleIdx/fftSize,sqrt(dftQ.^2+dftI.^2));
%% For comparison
sampledInput = sin(2*pi*f*Ts + sPhase);
Y = fft(sampledInput(1:1024),1024);
Pyy = Y.*conj(Y)/1024;
f = (0:1023);
figure;
plot(f,Pyy)
title('Power spectral density')
xlabel('Frequency (Hz)')
原因就在于两个不同的信号肯定会给你两个不同的频谱。查看下面的代码,你会发现你实际给出的dft算法的输入是sampledI+jsampledQ
。因此,您在这里所做的不仅仅是将原始信号分解为 In-phase and quadrature components
,而是在这里进行 Hilbert transform
——将 real
信号更改为 [=15] =]一个。
cI = cos(2*pi*w*sampleIdx/fftSize); %correlation cos
cQ = -sin(2*pi*w*sampleIdx/fftSize); %correlation sin
dftI(w+1) = sum(sampledI.*cI - sampledQ.*cQ); %
dftQ(w+1) = sum(sampledI.*cQ + sampledQ.*cI);
所以比较的 sampledInput
应该是 sampledInput = sI+1i*sQ;
.