Clear Filters
Clear Filters

Convolution problem

4 views (last 30 days)
yuvi
yuvi on 17 May 2011
hello to all
in this code we have an ofdm signal(60 MHZ),with length of 4096x1 (Pxx) and a power line model 0 to 80 Mhz with length of 1x161
the two parameters are shown in frequency domain
we have tried to get the plot after the ofdm signal is transmited threw the power line filter by convolution in time domain , or by multiplexing in frequensy domain
thank alot
the code :
clear all
%%%%%%%%%%%%%%%%%%%%%% O F D M S I G N A L %%%%%%%%%%%%%%%%%%%%%%
nFFTSize =64 ;
% for each symbol bits a1 to a52 are assigned to subcarrier
% index [-30 to -1 1 to 30]
subcarrierIndex = [-30:-1 1:30];
nBit = 2500;
ip1 = rand(1,nBit) >0.5;
ip2=rand(1,nBit) >0.5;
ip=[ip1;ip2] > 0.5; % generating 1's and 0's
nBitPerSymbol = 60;
nSymbol = ceil(nBit/nBitPerSymbol);
j=sqrt(-1)
% QBPSK modulation
% bit00 --> 1
% bit01 --> 1+j
% bit10 -->-1
% bit11 -->-1-j
t1 = ip1==1;
t2 = ip2==1;
ipMod(t1 & t2) = -(1+1i); %11
ipMod(~t1 & t2) = 1+1i; %01
ipMod(t1 & ~t2) = -1; %10
ipMod(~t1 & ~t2) = 1; %00
ipMod = [ipMod zeros(1,nBitPerSymbol*nSymbol-nBit)];
ipMod = reshape(ipMod,nSymbol,nBitPerSymbol);
figure(2)
plot(ipMod)
st = []; % empty vector
for ii = 1:nSymbol
inputiFFT = zeros(1,nFFTSize);
% assigning bits a1 to a52 to subcarriers [-30 to -1, 1 to 30]
inputiFFT(subcarrierIndex+nFFTSize/2+1) = ipMod(ii,:);
% shift subcarriers at indices [-30 to -1] to fft input indices [38 to 63]
inputiFFT = fftshift(inputiFFT);
outputiFFT = ifft(inputiFFT,nFFTSize);
% adding cyclic prefix of 16 samples
outputiFFT_with_CP = [outputiFFT(49:64) outputiFFT];
st = [st outputiFFT_with_CP];
end
close all
fsMHz = 60;
[Pxx,W] = pwelch(st,[],[],4096,20);
figure(3)
plot([-2048:2047]*fsMHz/4096,10*log10(fftshift(Pxx)));
xlabel('frequency, MHz')
ylabel('power spectral density')
title('Transmit spectrum OFDM (based on 802.11a)');
%%%%%%%%%%%%%%%%%%%%%%POWER LINE COMMUNICATION TRANSFER FUNCTION %%%%%%%%%%%%%%%%%%%%%%
w=0:0.5*10^6:80*10^6; % [Hz]
l=5 ; % [meter]
zl=50 ; % [ohm/meter]
zs=50 ; % [ohm/meter]
R=1.9884 ; % [ohm/meter]
G=0.01686*10^-9 ; % [mho/meter]
C=0.13394*10^-9 ; % [farad/meter]
L=362.81*10^-9 ; % [henrry/meter]
zc=sqrt((R+j.*w.*L)./(G+j.*w.*C)); % characteristic impedance [ohm]
gama=sqrt((R+j.*w.*L).*(G+j.*w.*C)); % propagation constant
a=cosh(gama.*l);
b=zc.*sinh(gama.*l);
c=(1./zc).*sinh(gama.*l);
d=cosh(gama.*l);
H= zl./((a.*zl)+b+(c.*zl.*zs)+(d.*zs)); %transfer function [Vl/Vs]
figure(1)
plot(w,20*log10(H)) % clean
xlabel('frequency [Hz]');
ylabel('Transfer function magnitude H(f) [dB]');
figure(5)
subplot(2,1,1)
plot(w,20*log10(H+0.001*[randn(1,length(H)) + j*randn(1,length(H))])) % plus random noise
xlabel('frequency [Hz] ');
ylabel('Transfer function + random noise [dB]');
subplot(2,1,2)
plot(w,20*log10(H+0.001*wgn(1,length(H),-10))) % plus white gaussian noise
xlabel('frequency [Hz] ');
ylabel('Transfer fuction + wgn [dB]');
y=conv(Pxx,H)
figure(4)
plot(20*log10(y))
  4 Comments
Arturo Moncada-Torres
Arturo Moncada-Torres on 17 May 2011
What do you want to fix from the axis?
yuvi
yuvi on 17 May 2011
in this command :
y=conv(Pxx,H)
figure(4)
plot(20*log10(y))
i want to get a plot with "x axis " of frequency
thanks

Sign in to comment.

Accepted Answer

Daniel Shub
Daniel Shub on 17 May 2011
Convolution of time domain signals is the same as multiplication of frequency domain signals. You cannot simply multiple Pxx and H since they have different sizes. You could take the ifft of both, conv and then take the fft of the result. Better would be to use W from pwelch instead of w in calculating H (you need to adjust for the sampling rate).
  3 Comments
Daniel Shub
Daniel Shub on 17 May 2011
Yes, you can just multiple in frequency, but NO you cannot just zero pad. You need to create H and Pxx with the same spectral spacing. Think about your W and w variables and what they represent.
To get the plot in the frequency domain try plot(W, 20*log10(y)).
On a side note, do you really want 20*log10 with Pxx (an estimate of power).
yuvi
yuvi on 17 May 2011
THANK'S DANIEL !

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!