plot_mimo.m

这是一个关于ofdm在802.11协议下的源码

% Plots the impulse responses, the PDPs, the correlation properties
% and the Doppler spectra of the MIMO channel stored in the global
% variable H produced by script example_MIMO.m. plot_MIMO.m also
% relies on other global variables of example_MIMO.m
%
% Revision history
%
% February 2004 - Bug fix related to the abscissa of the impulse
%                 response plot. Thanks to Jiun Siew, Bristol
%                 University (United Kingdom), for reporting it.
% November 2003 - Updated to rely on time-domain fading generation
% October 2003  - Updated in compliance with IEEE 802.11-03/161r2
%                 ( Bas Dijkstra )
%               - Modified bounds on Doppler spectrum to relate
%                 to Welch periodogram
% July 2003     - Modified to present IEEE 802-11 compliant results
% April 2003    - Modified to support more than 6 paths
%                 + addition of graphical check of Rayleigh
%                   distribution
%                 + addition of Welch PSD estimator with 50%
%                   overlap
%                 + bug fix related to the chip oversampling.
%                   Thanks to Terhi Rauhitiainen, Nokia Research
%                   Center, Helsinki (Finland) for reporting and
%                   fixing the bug.
% February 2002 - Creation
%
%
% STANDARD DISCLAIMER
%
% The Computer Science Institute of the University of Namur (hereafter
% "FUNDP-INFO") is furnishing this item "as is". FUNDP-INFO does not
% provide any warranty of the item whatsoever, whether express,
% implied, or statutory, including, but not limited to, any warranty
% of merchantability or fitness for a particular purpose or any
% warranty that the contents of the item will be error-free.
%
% In no respect shall FUNDP-INFO incur any liability for any damages,
% including, but not limited to, direct, indirect, special, or
% consequential damages arising out of, resulting from, or any way
% connected to the use of the item, whether or not based upon
% warranty, contract, tort, or otherwise; whether or not injury was
% sustained by persons or property or otherwise; and whether or not
% loss was sustained from, or arose out of, the results of, the
% item, or any services that may be provided by FUNDP-INFO.
%
% (c) Laurent Schumacher, FUNDP-INFO - February 2004

close all;

figure(1)
colour   = ['b','r','k','m','g','c'];
narrowH  = zeros(NumberOfTxAntennas*NumberOfRxAntennas,NumberOfHSamples);
abscissa = (0:1:NumberOfHSamples-1).*(DownSamplingFactor*SamplingTime_s);
summum   = find(PDP_linear(1,:) == max(PDP_linear(1,:)));
for ii = 1:NumberOfTxAntennas
    for jj = 1:NumberOfRxAntennas
        kk  = (ii-1)*NumberOfRxAntennas+jj;
        for ll = 1:NumberOfPaths
            tmp         = reshape(H(ii,jj,ll,:),1,NumberOfHSamples);
            narrowH(kk,:) = narrowH(kk,:)+tmp;
            subplot(NumberOfTxAntennas,NumberOfRxAntennas,kk),...
                semilogy(abscissa,abs(tmp),...
                colour(mod(ll,size(colour,2))+1)), ...
                hold on;
        end;
        grid;
        xlabel('Time [s]');
        title(['Tx#',num2str(ii),' - Rx#',num2str(jj)]);
        if (jj==1)
            subplot(NumberOfTxAntennas,NumberOfRxAntennas,kk),...
                ylabel('|H| [dB]');
        end;
    end;
end;

if (IEEE_802_11_Case == 'F')
    figure(2);
    for ii = 1:NumberOfTxAntennas
        for jj = 1:NumberOfRxAntennas
            kk  = (ii-1)*NumberOfRxAntennas+jj;
            subplot(NumberOfTxAntennas,NumberOfRxAntennas,kk),...
                semilogy(abscissa,abs(narrowH(kk,:))),...
                grid, xlabel('Time [s]'),...
                title(['Tx#',num2str(ii),' - Rx#',num2str(jj)]);
            if (jj==1)
                subplot(NumberOfTxAntennas,NumberOfRxAntennas,kk),...
                    ylabel('|\Sigma H| [dB]');
            end;
        end;
    end;
end;

% Check PDP

figure(3);
summum = find(PDP_linear(1,:) == max(PDP_linear(1,:)));
for ii = 1:NumberOfTxAntennas
    for jj = 1:NumberOfRxAntennas
        kk  = (ii-1)*NumberOfRxAntennas+jj;
        pdp = zeros(1,NumberOfPaths);
        for ll = 1:NumberOfPaths
            pdp(ll) = sum((abs(reshape(H(ii,jj,ll,:),1,NumberOfHSamples))).^2)./NumberOfHSamples;
        end;
        subplot(NumberOfTxAntennas,NumberOfRxAntennas,kk),...
            stem(1:NumberOfPaths,10*log10(PDP_linear(1,:)./PDP_linear(1,summum)),'r--'), hold on, ...
            stem(1:NumberOfPaths,10*log10(pdp./pdp(summum)),'b'), grid,...
            title(['Tx#',num2str(ii),' - Rx#',num2str(jj)]);
        axe    = axis;
        axe(1) = 1;
        axis(axe);
        if (jj==1)
            subplot(NumberOfTxAntennas,NumberOfRxAntennas,kk),...
                ylabel('Power [dB]');
        end;
        if (ii==NumberOfTxAntennas)
            subplot(NumberOfTxAntennas,NumberOfRxAntennas,kk),...
                xlabel('Tap index');
        end;
    end;
end;

% Check Rayleigh

figure(4);
h_ref = sqrt(0.5)*(randn(1,NumberOfHSamples)...
                   +randn(1,NumberOfHSamples)*j);
pdf_ref = hist(20*log10(abs(h_ref)),[-30:1:15]);
cdf_ref = cumsum(pdf_ref)./NumberOfHSamples;

% Warning disabled as log10(0) likely to occur
warning off;

for ii = 1:NumberOfTxAntennas
    for jj = 1:NumberOfRxAntennas
        for kk = 1:NumberOfPaths
            h_sim = reshape(H(ii,jj,kk,:),1,NumberOfHSamples);
            h_sim = h_sim./sqrt(mean(abs(h_sim).^2));
            pdf_sim = hist(20*log10(abs(h_sim)),[-30:1:15]);
            cdf_sim = cumsum(pdf_sim)./NumberOfHSamples;
            subplot(NumberOfTxAntennas,NumberOfRxAntennas,((ii-1)*NumberOfRxAntennas)+jj),...
                plot(-30:1:15,log10(cdf_sim),'b:'), hold on;
        end;
        subplot(NumberOfTxAntennas,NumberOfRxAntennas,((ii-1)*NumberOfRxAntennas)+jj),...
            plot(-30:1:15,log10(cdf_ref),'r'), grid,...
            title(['Tx#',num2str(ii),' - Rx#',num2str(jj)]);
        if (mod(jj,2) == 1)
            ylabel('log_{10} CDF');
        end;
        if (ii == NumberOfTxAntennas)
            xlabel('20 log_{10}(h) [dB]');
        end;
    end;
end;
warning on;

% Check correlation

if ((NumberOfTxAntennas*NumberOfRxAntennas) > 1)
    % No point in checking correlation of a SISO set-up
    Rcomplex = zeros(NumberOfPaths,NumberOfTxAntennas*NumberOfRxAntennas,NumberOfTxAntennas*NumberOfRxAntennas);
    for ii=1:NumberOfPaths
        for jj=1:NumberOfTxAntennas
            for kk=1:NumberOfRxAntennas
                ll = (jj-1)*NumberOfRxAntennas + kk;
                for mm=1:NumberOfTxAntennas
                    for nn=1:NumberOfRxAntennas
                        oo                 = (mm-1)*NumberOfRxAntennas + nn;
                        ha                 = reshape(H(jj,kk,ii,:),NumberOfHSamples,1);
                        hb                 = reshape(H(mm,nn,ii,:),NumberOfHSamples,1);
                        temp               = corrcoef(ha,hb);
                        Rcomplex(ii,ll,oo) = temp(1,2);
                    end;
                end;
            end;
        end;
    end;
    for ii = 1:NumberOfPaths
        figure(5 + floor((ii-1)/6));
        index_low  = (ii-1)*NumberOfTxAntennas*NumberOfRxAntennas+1;
        index_high = ii*NumberOfTxAntennas*NumberOfRxAntennas;
        for jj = 1:NumberOfTxAntennas*NumberOfRxAntennas
            for kk = 1:NumberOfTxAntennas*NumberOfRxAntennas
                temp(jj,kk) = abs(Rcomplex(ii,jj,kk)).^2;
            end
        end
        subplot(2,min(NumberOfPaths,6),1+mod((ii-1),6)),...
            mesh(1:NumberOfTxAntennas*NumberOfRxAntennas,1:NumberOfTxAntennas*NumberOfRxAntennas,temp),...
            title(['\langle h_{11}^{',num2str(ii),'}, h_{kl}^{',num2str(ii),'}\rangle']),...
        subplot(2,min(NumberOfPaths,6),min(NumberOfPaths,6)+1+mod((ii-1),6)),...
            plot(1:NumberOfTxAntennas*NumberOfRxAntennas, temp(1,:), 'b', ...
                 1:NumberOfTxAntennas*NumberOfRxAntennas, ...
                 abs(reshape(R(index_low:index_high,index_low), 1, ...
                 NumberOfTxAntennas*NumberOfRxAntennas)).^2, 'r--'),...
                 title(['\langle h_{11}^{',num2str(ii),'}, h_{kl}^{',num2str(ii),'}\rangle']),...
                 xlabel(['(k-1)*',num2str(NumberOfRxAntennas),' + l']);
        axe    = axis;
        axe(1) = 1;
        axe(3) = 0;
        axis(axe);
        grid;
        if (mod((ii-1),6) == 0)
            subplot(2,min(NumberOfPaths,6),1), zlabel('Correlation coefficient');
            subplot(2,min(NumberOfPaths,6),min(NumberOfPaths,6)+1), ylabel('Correlation coefficient');
        end;
    end;
end;

% Check Doppler spectrum

NumberOfDopplerSamples = 2^(floor(log2(NumberOfHSamples)));
switch FadingType
case {'bell_shape','bell_shape_spike'}
    normalised_frequency = (-NumberOfDopplerSamples/2+1:1:NumberOfDopplerSamples/2)...
                         * SamplingRate_Hz ...
                         / (DownSamplingFactor * NumberOfDopplerSamples * f_D_Hz);
otherwise
    normalised_frequency = (-NumberOfDopplerSamples/2+1:1:NumberOfDopplerSamples/2)...
                         * SamplingRate_Hz  ...
                         / (DownSamplingFactor * NumberOfDopplerSamples * Max_Doppler_shift);
end;
                   
% Parameters for Welch PSD estimation
% Length of a periodogram = 2^power of sample length - 4, at least 128
SectionLength = 2^(floor(log2(NumberOfDopplerSamples))-4);
if ((SectionLength < 128) & (NumberOfDopplerSamples > 128))
    SectionLength = 128;
end;
Downsampling = floor(NumberOfDopplerSamples/SectionLength);
% Number of periodograms with 50% overlap
NumberOfSections = floor(((2*NumberOfDopplerSamples)-SectionLength)/SectionLength);
% Time domain window
WelchWindow = hanning(SectionLength).';
% Normalisation factor as defined in Kunt M., "Traitement Numerique des Signaux", p. 267
NormalisationFactor = sum(WelchWindow.^2)/SectionLength;

for ii = 1:NumberOfTxAntennas
    for jj = 1:NumberOfRxAntennas
        for kk = 1:NumberOfPaths
            
            figure(5 + ceil(NumberOfPaths/6) + floor((kk-1)/6));
            tempH = reshape(H(ii, jj, kk, 1:NumberOfDopplerSamples), 1, NumberOfDopplerSamples);
            % Simple PSD estimator
            % (See Kunt M., "Traitement Numerique des Signaux", p. 256, Presses Polytechniques Romandes, 1984)
            spectrum = ((abs(fftshift(fft(tempH))).^2)./NumberOfDopplerSamples);
            
            % Welch PSD estimator with 50% overlap
            % (For Welch estimator, refer to Kunt M., "Traitement Numerique des Signaux", p. 256, Presses
            % Polytechniques Romandes, 1984; for overlap, refer to Carter G.C., Knapp C.H. and Nuttall A.H.,
            % "Estimation of the Magnitude-Squared Coherence Function via Overlapped Fast Fourier Transform
            % Processing", IEEE Transactions on Audio and Electroacoustics, vol. 21, n. 4, pp. 337-344,
            % August 1983)
            PeriodSpectrum = zeros(NumberOfSections,SectionLength);
            for ll = 1:NumberOfSections;
                SectionStart = 1+((ll-1)*SectionLength)/2;
                SectionEnd   = (ll+1)*SectionLength/2;
                PeriodSpectrum(ll,:) = ((abs(fftshift(fft(tempH(1,SectionStart:SectionEnd).*WelchWindow))).^2)./SectionLength);
            end;
            WelchSpectrum = sum(PeriodSpectrum,1)./(NumberOfSections*NormalisationFactor);
            
            switch FadingType
            case {'bell_shape','bell_shape_spike'}
                subplot(NumberOfTxAntennas*NumberOfRxAntennas,min(NumberOfPaths,6),...
                    (((ii-1)*NumberOfRxAntennas)+(jj-1))*min(NumberOfPaths,6)+mod((kk-1),6)+1),...
                    semilogy(normalised_frequency,spectrum,'b',...
                             normalised_frequency(1:Downsampling:NumberOfDopplerSamples),WelchSpectrum,'g');
            otherwise
                subplot(NumberOfTxAntennas*NumberOfRxAntennas,min(NumberOfPaths,6),...
                    (((ii-1)*NumberOfRxAntennas)+(jj-1))*min(NumberOfPaths,6)+mod((kk-1),6)+1),...
                    plot(normalised_frequency,spectrum,'b',...
                         normalised_frequency(1:Downsampling:NumberOfDopplerSamples),WelchSpectrum,'g');
            end;
            axe    = axis;
            switch FadingType
            case 'bell_shape_spike'
                axe(1) = max((-1)*ceil(1.25*(200/f_D_Hz)),axe(1));
                axe(2) = min(ceil(1.25*(200/f_D_Hz)),axe(2));
            otherwise
                axe(1) = max(-2,axe(1));
                axe(2) = min(2,axe(2));
            end;
            axe(3) = min(min(spectrum),min(WelchSpectrum));
            axe(4) = max(max(spectrum),max(WelchSpectrum));
            axis(axe);
            lower_bound = line(ones(1,2),[axe(3) axe(4)]);
            set(lower_bound,'Color',[1 0 0],'LineStyle',':');
            upper_bound = line(-1*ones(1,2),[axe(3) axe(4)]);
            set(upper_bound,'Color',[1 0 0],'LineStyle',':');
            switch FadingType
            case 'bell_shape_spike'
                if (kk==9)
                    spike_lower_bound = line((200/f_D_Hz)*ones(1,2),[axe(3) axe(4)]);
                    set(spike_lower_bound,'Color',[0 0 1],'LineStyle',':');
                    spike_upper_bound = line((-200/f_D_Hz)*ones(1,2),[axe(3) axe(4)]);
                    set(spike_upper_bound,'Color',[0 0 1],'LineStyle',':');
                end;
            end;
            switch FadingType
            case {'bell_shape','bell_shape_spike'}
                spectrum_upper_bound = line([axe(1) axe(2)],[max(WelchSpectrum) max(WelchSpectrum)]);
                set(spectrum_upper_bound,'Color',[0 1 0],'LineStyle',':');
                % Bound at 10 dB
                spectrum_lower_bound = line([axe(1) axe(2)],[.1*max(WelchSpectrum) .1*max(WelchSpectrum)]);
                set(spectrum_lower_bound,'Color',[0 1 0],'LineStyle',':');
            end;
            title(['Tap h_{',num2str(ii),num2str(jj),'}^{',num2str(kk),'}']);
        end;
    end;
end;