Convolution problem

2 Ansichten (letzte 30 Tage)
yuvi
yuvi am 17 Mai 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 Kommentare
Arturo Moncada-Torres
Arturo Moncada-Torres am 17 Mai 2011
What do you want to fix from the axis?
yuvi
yuvi am 17 Mai 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

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Daniel Shub
Daniel Shub am 17 Mai 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 Kommentare
Daniel Shub
Daniel Shub am 17 Mai 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 am 17 Mai 2011
THANK'S DANIEL !

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Community Treasure Hunt

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

Start Hunting!

Translated by