Reconstructing a Random Step Sequence using FFT

1 次查看(过去 30 天)
Dear All,
It may seem like a trivial question, but can't seem to figure it out.
I am using fft to reconstruct a randomly generated step sequence (0's and 1's) through Fourier series coefficients. All else is fine except my reconstructed signal is off by one sample, meaning it is shifted back by 1 sample. Can you please tell me what could I be doing wrong here?
Thank you for your time and help.
Regards, Adeel
P.S. you can skip to part 2 in the code where I am calculating the Fourier series coefficients.
%%Part 1
%%To obtain fourier series coefficients of a peroidic signal using fft
clear all; close all; clc;
dur = 0.5; % Signal duration (sec)
step = 1e-3; % step size
f = 50; % frequency of the signal (Hz)
t = 0:step:dur-step;
M = length(t);
y1 = zeros(1,M);
p = 1;
%%This for loop generates a random step sequence with different
% duty cycle each time it runs. It emulates the behavior of
% sinusoidal PWM for use in DC to AC converters.
for i = 1:M
n = randi([0,1],1); % generate a random 0 or 1
k = randi([10,50],1); % generate a random integer between 10 & 50
if n == 0; % if n is 0, then the k entries of y1 will be 0
for j = 1:k
y1(p) = 1;
p = p+1; % increment index of y1 vector
end
else % if n is 1, then the k entries of y1 will be 1
for j = 1:k
y1(p) = 0;
p = p+1; % increment index of y1
end
end
if p >=M
y1 = y1(1:M);
break;
end
end
figure(1)
plot(t,y1,'LineWidth',2)
axis([0 dur -0.1 1.1])
%%----------------------------------------------
%%Part 2
%%Calculate the Fourier series coefficients of y1
N = length(y1); % window length/data length
Y1 = fft(y1)/N; % take fft
N1 = ceil(N/2);
x = 1:1:N;
k = 0:N-1;
omega = k*2*pi*f; % fundamental frequencies and its harmonics
an = 2*real(Y1);
an(1) = an(1)/2;
bn = -2*imag(Y1);
fapp = an(1)*ones(size(y1));
L = N1;
for k = 1:L
fapp = fapp + an(k+1)*cos(2*pi*x*k/N) + bn(k+1)*sin(2*pi*x*k/N);
end
plot(t,y1,'-.b'); hold on
plot(t,fapp,'-r')

采纳的回答

Adeel
Adeel 2012-10-13
Resolved. I changed x from 0:1:N-1 instead of 1:1:N, and it fixed the issue.

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by