10 Referências bibliográficas - PUC-Rio · 10 Referências bibliográficas 1 FUHL, J. Smart antennas for second and third generation mobile communications systems. Viena - Áustria,

10 Referências bibliográficas

1 FUHL, J. Smart antennas for second and third generation mobile communications systems. Viena - Áustria, 1997. Ph.D dissertation - Technische Universität Wien.

2 REED, J.H. et al. An overview of the challenges and progress in meeting the E-911 requirement for location service. IEEE Communications Magazine, v.36, n.4, p. 30-37, Abr 1998.

3 CAFFERY, J.J.; STUBER, G.L. Overview of radiolocation in CDMA cellular systems. IEEE Communications Magazine, v.36, n.4, p. 38-45, Abr 1998.

4 DA ROCHA MARQUES, A. Localização geográfica de estações móveis em sistemas celulares. Rio de Janeiro, Set 2001. Dissertação de mestrado - Pontifícia Universidade Católica do Rio de Janeiro.

5 AHONEN, S.; ESKELINEN, P. Mobile terminal location for UMTS. IEEE Aerospace and Electronics Systems Magazine. v.18, n.2, p. 23-27, Fev 2003.

6 KUCHAR, A.; ROSSI, J.-P.; BONEK, E. Directional macro-cell channel characterization from urban measurements. IEEE Transactions on Antennas and Propagation, v.48, n.2, p. 137-146, Fev 2000.

7 THOMÄ, R. S. et al. Identification of time-variant directional mobile radio channels. In: IMTC, 16. 1999, Veneza - Itália. Proceedings...

8 MACÊDO, L. H. Sondagem em freqüência do canal indoor de faixa larga. Rio de Janeiro, Fev 2002. Dissertação de mestrado - Pontifícia Universidade Católica do Rio de Janeiro.

9 MACÊDO, L. H. et al. Mobile indoor wide-band 1.8 GHz sounding: measurement-based time dispersion analysis. In: IEEE Vehicular Technology Conference (VTC) spring 2002, 55. 2002, Birmingham – AL – EUA. Proceedings... v.1, p. 375-379.

10 VÁSQUEZ, E. J. A., Caracterização do canal móvel em faixa larga. Rio de Janeiro, Abr 2000. Tese de doutorado - Pontifícia Universidade Católica do Rio de Janeiro.

11 MACÊDO, J. F. Caracterização espaço-temporal do canal rádio móvel. Rio de Janeiro, Fev 2003. Dissertação de mestrado - Pontifícia Universidade Católica do Rio de Janeiro.

12 GODARA, L. C. Applications of antenna arrays to mobile communications, part I: performance improvement, feasibility, and system considerations. Proceedings of the IEEE, v.85, n.2, p. 1029-1060, Jul 1997.

13 ANDERSON, S. et al. An adaptive array for mobile communication systems. IEEE Transactions on Vehicular Technology, v.40, p. 230-236, 1991.

DBD

PUC-Rio - Certificação Digital Nº 0016247/CA

10 Referências bibliográficas 268

14 WINTERS, J. H.; SALZ, J.; GITLIN, R. D. The impact of antenna diversity on the capacity of wireless communication systems. IEEE Transactions on Communications, v.42, p. 1740-1751, 1994.

15 OHGANE, T. et al. BER performance of CMA adaptive array for high-speed GMSK mobile communication - A description of measurements in central Tokyo. IEEE Transactions on Vehicular technology, v.42, p. 484-490, 1993.

16 STOICA, P.; MOSES, R. Introduction to spectral analysis. Upper Saddle River, Prentice Hall PTR, 1997. 319p.

17 PISARENKO, V. F. The retrieval of harmonics from a covariance function. Geophysical Journal of the Royal Astronomy Society, v.33, p. 347-366, 1973.

18 SCHMIDT, R. O. Multiple emitter location and signal parameter estimation. IEEE Transactions on Antennas and Propagation, v.34, p. 276-280, 1986.

19 ROY, R.; KAILATH, T. ESPRIT – Estimation of Signal Parameters via Rotational Invariance techniques. IEEE Transactions on Acoustics, Speech, and Signal Processing, v.37, n.7, p. 984-995, Jul 1989.

20 JOHNSON, D. H. The application of spectral estimation methods to bearing estimation problems. Proceedings of the IEEE, v.70, p. 1018-1028, 1982.

21 CAPON, J. High-resolution frequency-wave number spectrum analysis. Proceedings of the IEEE, v.57, p. 1408-1418, 1969.

22 LIBERTI JR, J. C.; RAPPAPORT, T. S. Smart antennas for wireless communications: IS-95 and third generation CDMA applications. Upper Saddle River: Prentice Hall PTR, 1999.

23 GODARA, L. C. Applications of antenna arrays to mobile communications, part II: beam-forming and direction-of-arrival considerations. Proceedings of the IEEE, v.85, n.8, p. 1193-1245, Ago 1997.

24 PAULRAJ, A.; PAPADIAS, C. B. Array processing for mobile communications. In:___Handbook on signal processing. CRC Press, 1997.

25 ROY, R.; KAILATH, T. ESPRIT – Estimation of Signal Parameters via Rotational Invariance Techniques. IEEE Tranactions on Acoustics, Speech and Signal Processing, v.37, n.7, p. 984-995, Jul 1989.

26 RAPPAPORT, T.S. Wireless communications: principles and practice. New Jersey: Prentice-Hall, 1996. 641p.

27 PARSONS, D. The mobile radio propagation channel. John Wiley & Sons: New York-Toronto, 1992, ch. 8, pp. 212-254.

28 FUHL, J.; ROSSI, J.-P.; BONEK, E. High-resolution 3-D direction-of-arrival determination for urban mobile radio. IEEE Transations on Antennas and Propagation, v.45, n.4, p. 672-682, Abr 1997.

29 ROSSI, J. et al. Theory and measurement of the angle of arrival and time delay of UHF radiowaves using a ring array. IEEE Transations on Antennas and Propagation, v.45, n.5, p. 876-884, Mai 1997.

DBD



30 HU, Z. Evaluation of joint aoa and doa estimation algorithms using the antenna array systems. Blacksburg – Virginia – EUA, Dez 1998. MSc. Thesis - Virginia Polytechnic Institute and State University.

31 QIAN, S.; CHEN, D. Joint time-frequency analysis: methods and applications. Upper Saddle River – NJ – EUA: Prentice Hall PTR, 1996.

32 WICKERHAUSER, M. V. Adapted wavelet analysis: from theory to software. Piscataway – NJ – EUA: IEEE Press, 1994.

33 MISITI, M.; OPPENHEIM, G.; POGGI, J.-M. Matlab 6.0 wavelet toolbox help. The Mathworks Incorporation, Set 2000.

34 XU, W.; STEWART, W. K. Multiresolution-signal direction-of-arrival estimation: a wavelets approach. IEEE Transactions on Signal Processing Letters, v.7, n.3, p. 66-68, Mar 2000.

35 MANTIS, S.; HIPPENSTIEL, R. Time difference of arrival estimation of unequal SNR denoised communication signals. In: Asilomar Conference on Signals, Systems and Computers, 34. 2000. Conference Record... v.1, p. 629-633.

36 CHING, P.C.; SO, H.C.; WU, S.Q. On wavelet denoising and its applications to time delay estimation. IEEE Transactions on Signal Processing, v.47, n.10, p. 2879-2882, Out 1999.

37 NIU, X.X.; CHING, P.C.; CHAN, Y.T. Wavelet based approach for joint time delay and Doppler stretch measurements. IEEE Transactions on Aerospace and Electronic Systems, v.35, n.3, p. 1111-1119, Jul 1999.

38 HIPPENSTIEL, R.; AKTAS, U. Difference of arrival estimation using wavelet transforms. In: Midwest Symposium on Circuits and Systems, 42. 2000. Proceedings... v.2 , p. 1082-1085.

39 RODRIGUEZ, C.; BRAVO, A.M.; BHARGAVA, V.K. Detection and delay estimation of a new user entering in a CDMA system. In: IEEE Vehicular Technology Conference (VTC), 48. 1998. Proceedings... v.3, p. 1730-1734.

40 WU, S.Q.; SO, H.C.; CHING, P.C. Improvement of TDOA measurement using wavelet denoising with a novel thresholding technique. IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP). 1997. Proceedings... v.1, p. 539-542.

41 QUI, A.; BOULINGUEZ, D. Multipath channel identification with wavelet packets. IEEE Journal of Oceanic Engineering, v.22, n.2, p. 342-346, Abr 1997.

42 FARGUES, M. P. et al. Wavelet-based denoising: comparisons between orthogonal and non-orthogonal decompositions. In: Midwest Symposium on Circuits and Systems, 40. 1997. Proceedings... v.2, p. 929-932.

43 DONOHO, D. L.; JOHNSTONE, I. M. Ideal spatial adaptation via wavelet shrinkage. Biometrika, v.81, p. 425-455, Set 1994.

44 TASWELL, C. The what, how, and why of wavelet shrinkage denoising. Computing in Science & Engineering, v.2, n.3, p. 12-19, Mai-Jun 2000.

45 DONOHO, D. L. De-noising by soft-thresholding. IEEE Transactions on Information Theory, v.41, n.3, p. 613-627, Mai 1995.

DBD



46 HUETHER, B. M.; GUSTAFSON, S. C.; BROUSSARD, R. P. Wavelet preprocessing for high range resolution radar classification. IEEE Transactions on Aerospace and Electronic Systems, v.37, n.4, p. 1321-1332, Out 2001.

47 SOUSA, E. S.; JOVANOVIC, V. M.; DAIGNEAULT, C. Delay spread measurements for the digital cellular in Toronto. IEEE Transactions on Vehicular Technology, v.43, p. 837-847, Nov 1994.

48 HOWARD, S. T.; PAHLAVAN, K. Measurement and analysis of the indoor radio channel in frequency domain. IEEE Transactions of Instrumentation and Measurement, v.39, p. 751-755, Out 1990.

49 HARRIS, F. J. On the use of windows for harmonic analysis with the discrete Fourier transform. IEEE Proceedings, v.66, pp. 51-83, 1978.

50 DAUBECHIES, I. The wavelet transform, time-frequency localization and signal analysis. IEEE Transactions on Information Theory, v.36, n.5, p. 961-1005, Set 1990.

51 FUCHS, J.-J. Estimating the number of sinusoids in additive white noise. IEEE Transactions on Acoustics, Speech, and Signal Processing, v.36, n.12, p. 1846-1853, Dez 1988.

52 FUCHS, J.-J. Estimating the number of signals in the presence of unknown correlated sensor noise. IEEE Transactions on Signal Processing, v.40, n.5, p. 1053-1061, Mai 1992.

53 BALANIS, C. A. Advanced engineering electromagnetics. New York: John Wiley & Sons, 1989. 981p.

54 DIAS, M. H. C. Análise crítica da propagação em microcélulas. Rio de Janeiro, Jan 1998. Tese de mestrado – Instituto Militar de Engenharia.

55 DIAS, M. H. C.; ASSIS, M. S. Some remarks on urban and suburban microcell propagation. In: URSI Commission F International Open Symposium on Wave Propagation and Remote Sensing. 1998, Aveiro - Portugal. Proceedings... p. 224–227.

56 DIAS, M. H. C.; RAMOS, G. L.; SIQUEIRA, G. L. Ray-tracing analysis of 3.5 GHz propagation at a typical urban environment. SBMO/IEEE MTT-S International Microwave and Optoelectronics Conference (IMOC). 2001, Belém. Proceedings... v.1, p. 203–207.

57 MENENDEZ, R. C.; LEE, S. W. On the role of the geometrical optics field in aperture diffraction. IEEE Transactions on Antennas and Propagation, v. AP-25, n.5, p. 688-695, Set 1977.

58 TINGLEY, R. D.; PAHLAVAN, K. Space-time measurement of indoor radio propagation. IEEE Transactions on Instrumentation and Measurements, v.50, p. 22-31, Fev 2001.

59 DAL BELLO, J. C. R. Caracterização da influência da vegetação nos sistemas de comunicações móveis celulares em áreas urbanas. Rio de Janeiro, Fev 1998. Tese de doutorado - Pontifícia Universidade Católica do Rio de Janeiro.

DBD



60 SATISH, A. et al. Maximum likelihood estimation and Cramer-Rao bounds for direction of arrival parameters of a large sensor array. IEEE Transactions on Antennas and Propagation, v.44, n.4, p. 478-490, Abr 1996.

61 ZELENOVSKY, R. Emprego de arranjo de antenas na recuperação de dados digitais em ambiente CDMA. Rio de Janeiro, Jul 2001. Tese de doutorado - Pontifícia Universidade Católica do Rio de Janeiro.

62 ZOLTOWSKI, M. D.; HAARDT, M.; MATHEWS, C. P. Closed-form 2-D angle estimation with rectangular arrays in element space or beamspace via unitary ESPRIT. IEEE Transactions on Signal Processing, v.44, n.2, p. 316-328, Fev 1996.

DBD


Apêndice A Algoritmos utilizados na aplicação da supressão de ruído a PDPs

Foram utilizados diversos algoritmos ao longo do trabalho. Eventualmente,

pequenas modificações eram realizadas para modificar os parâmetros de entrada

ou as formas de apresentação e armazenamento dos dados. Com isso, apenas

alguns dos algoritmos serão apresentados integralmente aqui, sem perda de

representatividade para o conjunto.

O algoritmo a seguir realiza a supressão de ruído no domínio da freqüência

em escala linear, para uma das 13 rotas de sondagem de [8], escolhida como

parâmetro de entrada. É importante destacar que os algoritmos usados em [8]

foram usados como referência para o desenvolvimento das rotinas de estimação

dos PDPs nesta tese. Denoising 1d indoor na freq. % versao de 12/02/03 % ------------------------------------------------------------------------------------- clear; clc; %limpa tudo addpath ..\origem\; addpath ..\resultados\; addpath ..\funcoes\; disp ('Denoising 1D no dominio da freqüência, antes do janelamento - esc lin'); hs=input('s ou h: ','s'); level=input('nr de decomposicoes (max 5): '); arq=input('rota: ','s'); rota=['c1','c2','c3','l1','g1','g2','g3','g4','g5','g6','g7','g8','g9']; drota=[38,8,35,140,51,18,47,40,20,54,20,41,73]; for v=1:13 if rota((2*v-1):(2*v))==arq distancia=drota(v);disp(distancia); end end span= 200; disp('Carregando arquivos...'); arqm= sprintf('%sM.txt',arq); magn=load(arqm); arqf=sprintf('%sF.txt',arq); fase=load(arqf); clear arqm arqf; arqt=sprintf('%sT.txt',arq); tempo=load(arqt); clear arqm arqf arqt; disp('Arquivos carregados.'); [varreduras, pontos]=size(magn); %detecta tamanho %------------- %DEFINICOES DO DENOISING - MEXER SO AQUI!!! selec='sqtwolog'; %sqtwolog ou rigrsure scalam='mln'; scalas='sln'; %one ou sln ou mln wave='sym8'; comp=sprintf('df%s%sL%d',hs,wave,level); disp(comp); %-------------------------------------------------------------------------------------- dresolu=distancia/(varreduras-1); dist=0:dresolu:distancia; %--------------------------------------------------------------------------------------

DBD


Apêndice A Algoritmos utilizados na aplicação da supressão de ruído a PDPs 273

% Aproxima 'pontos' para adequar o dominio de retardos ao denoising, truncando o final % Para isso, o nr de amostras deve ser multiplo de 2^level %-------------- pontos=floor(pontos/(2^level))*(2^level); magn=magn(:,1:pontos); fase=fase(:,1:pontos); %-------------------------------------------------------------------------------------- % denoising disp('Realizando denoising sobre a amplitude do sinal...'); freq=(1:pontos)*span/pontos+1700; magn=10.^(magn/20); %converte p esc lin - campo for j=1:varreduras magnds(j,:)=wden(magn(j,:),selec,hs,scalas,level,wave); magndm(j,:)=wden(magn(j,:),selec,hs,scalam,level,wave); zd(j,:)=magn(j,:).*cos(fase(j,:)*pi/180)+i*magn(j,:).*sin(fase(j,:)*pi/180); zds(j,:)=magnds(j,:).*cos(fase(j,:)*pi/180)+i*magnds(j,:).*sin(fase(j,:)*pi/180); zdm(j,:)=magndm(j,:).*cos(fase(j,:)*pi/180)+i*magndm(j,:).*sin(fase(j,:)*pi/180); plot (freq,magn(j,:),'c');hold on; plot(freq,magnds(j,:),':b','LineWidth',2); plot(freq,magndm(j,:),'r','LineWidth',2); %title(sprintf('route: %s -- scheme: %slin -- %4.2f out of %4.2f m ran',arq,comp,dist(j),dist(varreduras))); title(sprintf('rota: %s -- esquema: %slin -- %4.2f de %4.2f m percorridos',arq,comp,dist(j),dist(varreduras))); ylabel('E [V/m]'); xlabel ('f [MHz]'); legend('STR','SLN', 'MLN'); pause; hold off; end disp('Denoising realizado.'); clear magn fase magnds magnds; %-------------------------------------------------------------------------------------- % Janelamento w=BHarrisW(pontos); %funcao janela for j=1:varreduras %aplica a funcao janela zdj(j,:)=zd(j,:).*w; zdsj(j,:)=zds(j,:).*w; zdmj(j,:)=zdm(j,:).*w; end clear w zd zds zdm; %------------------------------ % Estimativa do PDP por correlograma disp('Estimando correlogramas...'); for j=1:varreduras aux2=ccorrelogramse(zdj(j,:), pontos); aux2=fliplr(aux2'); result(j,:)=20*log10(abs(aux2)); clear aux2; aux2=ccorrelogramse(zdsj(j,:), pontos); aux2=fliplr(aux2'); resultds(j,:)=20*log10(abs(aux2)); clear aux2; aux2=ccorrelogramse(zdmj(j,:), pontos); aux2=fliplr(aux2'); resultdm(j,:)=20*log10(abs(aux2)); clear aux2; end disp('Correlogramas estimados.'); %------------------------------------------------------------------------------------- retardo_max=10^3/(span/(pontos-1)); resolucao=10^3/span; %nano segundos retardo=0:resolucao:retardo_max; % ------------------------------------------------------------------------------------- limiret=2000; pt_limiret=round(limiret/resolucao); % retardo=retardo(1:pt_limiret); validos=result(:,1:pt_limiret); validosds=resultds(:,1:pt_limiret); validosdm=resultdm(:,1:pt_limiret); clear result resultds resultdm ; % ------------------------------------------------------------------------------------- disp('Calculando a media entre perfis...'); media=mean(10.^(validos/10)); media=10*log10(media); mediads=mean(10.^(validosds/10)); mediads=10*log10(mediads); mediadm=mean(10.^(validosdm/10)); mediadm=10*log10(mediadm); figure; plot(retardo,media,'c'); hold on; plot(retardo,mediads,':b','LineWidth',2); plot(retardo,mediadm,'r','LineWidth',2); legend('STR','SLN','MLN'); grid on; %title(sprintf('route: %s -- scheme: %s - linear scale',arq,comp)); title(sprintf('rota: %s -- esquema: %s - linear',arq,comp)); ylabel('P [dB]'); xlabel ('\tau [ns]'); hold off; arqfigm=sprintf('%s%slinmed',arq,comp); %grava fig media em jpg clear media mediads mediadm;

DBD



O algoritmo da função janela de Blackman-Harris (BHarrisW) foi extraído

de [8], e é dado por: function [Vetor]=BHarrisW (Nptos) % Esta Função implementa uma janela de no mínimo 3 termos % de Blackman-Harris. % Nptos é um parâmetro de entrada, indica o tamanho da janela, % e deve ser impar. % Vetor é a saída da Função % sendo que os coeficientes são a0=0.42323; a1=0.49755; a2=0.07922; for n=0:(Nptos-1) Vetor(n+1)=a0-a1*cos(2*pi*n/Nptos)+a2*cos(2*pi*2*n/Nptos); end

Uma diferença entre os algoritmos desta tese e os de [8] foi o uso da

estimação por correlograma ao invés da transformada inversa de Fourier discreta.

As funções necessárias foram extraídas de [16], mas precisaram ser adaptadas

para fornecer a estimativa complexa ao invés de seu valor absoluto. Na verdade

esta modificação foi mais voltada para os algoritmos de estimação do espectro

espacial-temporal; para a estimação do espectro temporal, bastaria o módulo. Os

algoritmos correspondentes às funções chamadas (cperiodogramse e

ccorrelogramse) são apresentados a seguir: function phi=cperiodogramse(y,v,L) % phi=periodogramse(y,v,L) % y -> vetor de dados % v -> vetor janela % L -> nr de amostras da resposta espectral % phi <- estimativas espectrais em L freqüências w=0, 2*pi/L, ..., 2*pi(L-1)/L % Modificado para fornecer a estimativa complexa ao inves do modulo simplesmente (2/8/2002 - MHCDias) % Para isto, foi retirado o operador abs(...).^2, que deve ser usado no algoritmo que chamar esta funcao % checa a largura do vetor janela M=length(v); N=length(y); if (M>N) error('A largura da janela é maior que a do vetor de dados'); return elseif (M<N), fprintf('ALERTA: A largura da janela é menor que a largura\n') fprintf(' do vetor de dados, que serah truncado\n') fprintf(' para a largura da janela\n') end y=y(:); phi=(fft((y(1:M).*v(:)),L))/M;

function phi=ccorrelogramse(y,L) % phi=correlogramse(y,L) % y -> vetor de dados % L -> nr de amostras da resposta espectral % phi <- estimativas espectrais em L freqüências w=0, 2*pi/L, ..., 2*pi(L-1)/L % Modificado para forncer a resposta complexa ao inves do modulo (2/8/2002 - MHCDias) y=y(:); phi=cperiodogramse(y,ones(size(y)),L);

DBD



Para exemplificar a supressão 2D no domínio do retardo, foi escolhido o

algoritmo a seguir, que atua de uma só vez sobre todas os dados coletados em [8]. % Denoising 2d indoor em PDP % versao de 12/02/03 % ------------------------------------------------------------------------------------- clear; clc; %limpa tudo addpath ..\origem\; addpath ..\result_den2\; addpath ..\funcoes\; %scala=input('sln ou mln ou one: ','s'); %hs=input('s ou h: ','s'); level=input('nr de decomposicoes (max 5): '); rota=['c1','c2','c3','l1','g1','g2','g3','g4','g5','g6','g7','g8','g9']; drota=[38,8,35,140,51,18,47,40,20,54,20,41,73]; for v=1:13 disp ('Denoising 2D sobre os PDPs - esc linear'); arq = rota((2*v-1):(2*v));disp(arq);%input('Nome da Medida : ','s'); span= 200; distancia=drota(v);disp(distancia);%input('distancia maxima [m]: '); disp('Carregando arquivos.. ..'); arqm= sprintf('%sM.txt',arq); magn=load(arqm); arqf=sprintf('%sF.txt',arq); fase=load(arqf); clear arqm arqf; arqt=sprintf('%sT.txt',arq); tempo=load(arqt); clear arqm arqf arqt; disp('Arquivos carregados.. ..'); [varreduras, pontos]=size(magn); %detecta tamanho pontos=800; %trunca o ´ultimo %------------- %DEFINICOES DO DENOISING - MEXER SO AQUI!!! %level=5; %nr de niveis de decomposicao para o denoising selec='sqtwolog'; %sqtwolog ou rigrsure scalam='mln'; scalas='sln'; %one ou sln ou mln wave='sym8'; hs='s'; %s ou h comp=sprintf('dpc%s%sL%d',hs,wave,level); disp(comp); %-------------------------------------------------------------------------------------- dresolu=distancia/(varreduras-1); dist=0:dresolu:distancia; %-------------------------------------------------------------------------------------- magn=10.^(magn/20); %corrigido para/20 em 19/08/2002, pois S21 do an rede eh relacao de POTENCIA!!! disp('Compondo o sinal complexo para janelamento e fft'); for j=1:varreduras %compoe o numero complexo no dominio original absoluto for p=1:pontos z(j,p)=magn(j,p)*cos(fase(j,p)*pi/180)+i*magn(j,p)*sin(fase(j,p)*pi/180); end end clear magn fase; %-------------------------------------------------------------------------------------- % Janelamento w=BHarrisW(pontos); %funcao janela for j=1:varreduras %aplica a funcao janela zj(j,:)=z(j,:).*w; end clear w z; %------------------------------ % Estimativa do PDP por correlograma disp('Estimando correlogramas...'); for j=1:varreduras aux2=ccorrelogramse(zj(j,:),pontos); aux2=fliplr(aux2'); aux3=angle(aux2); fase(j,:)=aux3(1,:); result(j,:)=20*log10(abs(aux2)); resultrs(j,:)=wden(abs(aux2).*cos(fase(j,:)),selec,hs,scalas,level,wave); resultis(j,:)=wden(abs(aux2).*sin(fase(j,:)),selec,hs,scalas,level,wave); resultds(j,:)=20*log10(abs(resultrs(j,:)+i*resultis(j,:))); resultrm(j,:)=wden(abs(aux2).*cos(fase(j,:)),selec,hs,scalam,level,wave); resultim(j,:)=wden(abs(aux2).*sin(fase(j,:)),selec,hs,scalam,level,wave); resultdm(j,:)=20*log10(abs(resultrm(j,:)+i*resultim(j,:))); clear aux2 aux3;

DBD



end clear fase resultrs resultis resultrm resultim; disp('Denoising realizado.'); retardo_max=10^3/(span/(pontos-1)); resolucao=10^3/span; %nano segundos retardo=0:resolucao:retardo_max; limiret=2000; pt_limiret=round(limiret/resolucao); retardo=retardo(1:pt_limiret); validos=result(:,1:pt_limiret); validosds=resultds(:,1:pt_limiret); validosdm=resultdm(:,1:pt_limiret); clear result resultds resultdm; % ------------------------------------------------------------------------------------- disp('Calculando a media entre perfis...'); media=mean(10.^(validos/10)); media=10*log10(media); mediads=mean(10.^(validosds/10)); mediads=10*log10(mediads); mediadm=mean(10.^(validosdm/10)); mediadm=10*log10(mediadm); figure; plot(retardo,media,'c'); hold on; plot(retardo,mediads,':b','LineWidth',2); plot(retardo,mediadm,'r','LineWidth',2); legend('STR','SLN','MLN'); grid on; %title(sprintf('route: %s -- scheme: %s - linear scale',arq,comp)); title(sprintf('rota: %s -- esquema: %s - linear',arq,comp)); ylabel('P [dB]'); xlabel ('\tau [ns]'); hold off; arqfigm=sprintf('%s%slinmed',arq,comp); %grava fig media em jpg cd ..\result_den2\figuras; print('-dbmp',arqfigm); cd ..\..\alg_den2; clear media mediads mediadm; disp('Salvando dados para analise .. ..'); destino=sprintf('..%sresult_den2%s%s%slin_stm','\','\',arq,comp); f1='dist'; f2='resolucao'; f4='validos'; f5='validosds'; f6='validosdm'; save(destino, f1, f2, f4,f5,f6); clear destino dist resolucao validos validosds validosdm; end

O algoritmo a seguir foi utilizado para calcular os momentos estatísticos

associados à faixa dinâmica e à relação sinal/ruído das medidas, para a supressão

de ruído 2D. clc; clear; addpath ..\origem\; addpath ..\result_den2\; addpath ..\funcoes\; addpath ..\alg_den2\; i1=input('dfc ou dpc: ','s'); i2='s'; i4='sym8'; %input('wave: ','s'); i5=input('decomposition levels: '); i6=input ('lin ou db: ','s'); disp('Carregando '); rota=['c1','c2','c3','l1','g1','g2','g3','g4','g5','g6','g7','g8','g9']; arq3=[i1 i2 i4 'L' int2str(i5) i6]; disp(arq3); for v=1:13 arq1 = rota((2*v-1):(2*v)); disp(arq1); arqm2= sprintf('%s%s_stm.mat',arq1,arq3); load(arqm2); pdp=validos; %pdp: pdp str pdpds=validosds; pdpdm=validosdm; clear validos validosds validosdm; dimr=size(pdp,2);dimd=size(pdp,1); retardo=[0:5:((dimr-1)*resolucao)]; pdp=10.^(pdp/10); pdpds=10.^(pdpds/10); pdpdm=10.^(pdpdm/10); disp('carregados'); %-------------------------------------------------------------------------- %Definicao do ponto inicial para calculo do piso de ruido mediano - rtiruido %Valores determinados por observacao dos pefis medios STR por percurso %MODIFICADO EM 18/02/2003 rtiruido=400; %g3 if (v==5)|(v==6)|(v==8)|(v==11) rtiruido=700; end if (v==12)|(v==9)|(v==13)

DBD



rtiruido=1000; end if (v==10)|(v==4)|(v==1)|(v==2)|(v==3) rtiruido=1200; end rtiruido=rtiruido/resolucao; %rtiruido eh o indice associado ao retardo %-------------------------------------------------------------------------- %estatisticas %DR e SNR %----pdp----------- dr=max(pdp')./mean(pdp(:,rtiruido:dimr)'); dr=10*log10(mean(dr)); snr=max(pdp')./(std(pdp(:,rtiruido:dimr),1,2))'; snr=10*log10(mean(snr)); %----pdpds--------- drs=max(pdpds')./mean(pdpds(:,rtiruido:dimr)'); drs=10*log10(mean(drs)); snrs=max(pdpds')./(std(pdpds(:,rtiruido:dimr),1,2))'; snrs=10*log10(mean(snrs)); %----pdpdm--------- drm=max(pdpdm')./mean(pdpdm(:,rtiruido:dimr)'); drm=10*log10(mean(drm)); snrm=max(pdpdm')./(std(pdpdm(:,rtiruido:dimr),1,2))'; snrm=10*log10(mean(snrm)); %----pdp - pdpds--- disp('STR x SLN'); disp('Valores medios para as faixas dinamicas'); textopisodr=sprintf('STR: %4.2f dB --- SLN: %4.2f dB --- dif: %4.2f',dr, drs,drs-dr); disp(textopisodr); disp('Valores medios para as SNR'); textopisosnr=sprintf('STR: %4.2f dB --- SLN: %4.2f dB --- dif: %4.2f',snr, snrs, snrs-snr); disp(textopisosnr); %----pdp - pdpdm--- disp('STR x MLN'); disp('Valores medios para as faixas dinamicas'); textopisodr=sprintf('STR: %4.2f dB --- MLN: %4.2f dB --- dif: %4.2f',dr, drm,drm-dr); disp(textopisodr); disp('Valores medios para as SNR'); textopisosnr=sprintf('STR: %4.2f dB --- MLN: %4.2f dB --- dif: %4.2f',snr, snrm, snrm-snr); disp(textopisosnr); clear pdp pdpds pdpdm dr drs drm snr snrs snrm arqm2; end;

Os dados coletados em [10] só foram processados pela técnica 1D, no

domínio do retardo, e em escala linear. Os algoritmos a seguir foram utilizados

para aplicar a supressão de ruído aos dados coletados pelas sondas MF e STDCC,

respectivamente. Em ambos, o usuário é solicitado a avaliar se cada resposta é

válida ou não. Apenas as respostas aprovadas são armazenadas para análise

posterior. O algoritmo para os dados da sonda STDCC é mostrado apenas

parcialmente, pois o segmento de armazenamento das respostas é análogo ao da

sonda MF.

%Captura dados gerados pelo programa ANALISE da tese do Eduardo clear; clc; cwd = pwd; cd(tempdir); pack; cd(cwd); %limpa tudo addpath ..\origem\MF\cetuc\; addpath ..\origem\MF\copa\; disp ('Programa de Captura dos resultados do programa ANALISE da tese do Eduardo'); arq = input('Nome do Arquivo : ','s'); disp('Carregando arquivo.. ..'); arqres= sprintf('%s.res',arq); fid=fopen(arqres); c=1; N=4096; %Nr de pontos por perfil try

DBD



while (ftell(fid)~=-1) lixo=fscanf(fid,'%s',[1 14]); %Le primeiras linhas de string ateh o valor inicial de retardo disp (sprintf('Varredura %u',c)); T(c)=fscanf(fid,'%lg',[1 1]); %Le valor inicial de retardo lixo2=fscanf(fid,'%s',[1 1]); %Le 'us' for i=1:N pdp(c,i)=fscanf(fid,'%lg',[1 1]); end c=c+1; clc; end catch disp('Fim do arquivo'); c=c-1; end fclose(fid); %-------------------------------------------------------------------------------------- % Define vetor de retardo: resposta com 51,1 us em 2555 ptos N2=2555; retardo=((1:N)-1)*(51.1/N2); %Define resposta em escala linear pdp=10.^(pdp/20); pdp=pdp(:,1:N2); retardo=retardo(1:N2); %------------------------------------- % Define parametros do denoising sh='s'; onda='sym8'; esqlimiar='sqtwolog'; L=7; parden=sprintf('%sL%d',sh,L); %------------------------------------- % Calcula resposta denoised na escala linear disp('Realizando denoising sobre os perfis em escala linear...'); figure; jc=1;nc=0; for j=1:c %denoising pdpds(j,:)=abs(wden(pdp(j,:),esqlimiar,sh,'sln',L,onda)); pdpdm(j,:)=abs(wden(pdp(j,:),esqlimiar,sh,'mln',L,onda)); minimo=20*log10(min(pdp(j,:))); if (minimo>-10) eixomin=-10; else eixomin=10*floor(minimo/10); end %Calcula SNR e DR snr(j)=20*log10(max(pdp(j,:))/std(pdp(j,:),1,2)); snrds(j)=20*log10(max(pdpds(j,:))/std(pdpds(j,:),1,2)); snrdm(j)=20*log10(max(pdpdm(j,:))/std(pdpdm(j,:),1,2)); dr(j)=20*log10(max(pdp(j,:))/mean(pdp(j,:))); drds(j)=20*log10(max(pdpds(j,:))/mean(pdpds(j,:))); drdm(j)=20*log10(max(pdpdm(j,:))/mean(pdpdm(j,:))); %Poe em escala log e plota pdpds(j,:)=20*log10(pdpds(j,:)); pdpdm(j,:)=20*log10(pdpdm(j,:)); pdp(j,:)=20*log10(pdp(j,:)); plot (retardo+T(j),pdp(j,:),'c');hold on; plot(retardo+T(j),pdpds(j,:),':b','LineWidth',2); plot(retardo+T(j),pdpdm(j,:),'r','LineWidth',2); title(sprintf('Perfil %u de %u da rota %s',j,c,arq)); xlabel('\tau [us]'); ylabel('P [dB]');grid on; axis([T(j) T(j)+50 eixomin 0]); legend('STR','SLN','MLN'); hold off; ask=input('valida? ','s'); if (ask~='n')&(ask~='N') pdpc(jc,:)=pdp(j,:);pdpdsc(jc,:)=pdpds(j,:);pdpdmc(jc,:)=pdpdm(j,:); snrc(jc)=snr(j);snrdsc(jc)=snrds(j);snrdmc(jc)=snrdm(j); drc(jc)=dr(j);drdsc(jc)=drds(j);drdmc(jc)=drdm(j); jc=jc+1;nc=jc-1; end end disp([int2str(c-nc) ' respostas cortadas']); clear pdp pdpds pdpdm snr snrds snrdm dr drds drdm; pdp=pdpc; pdpds=pdpdsc; pdpdm=pdpdmc; snr=snrc; dr=drc; snrds=snrdsc; drds=drdsc; snrdm=snrdmc; drdm=drdmc; clear pdpc pdpdsc pdpdmc snrc snrdsc snrdmc drc drdsc drdmc; if nc>1 %Calcula medias e grava vetores mdr=20*log10(mean(10.^(dr/20))); msnr=20*log10(mean(10.^(snr/20))); mdrds=20*log10(mean(10.^(drds/20))); msnrds=20*log10(mean(10.^(snrds/20))); mdrdm=20*log10(mean(10.^(drdm/20))); msnrdm=20*log10(mean(10.^(snrdm/20))); arqsnr=sprintf('%s%s-snr.xls',arq,parden); cd ..\results2\MF\planilhas; fid2 = fopen(arqsnr,'w'); fprintf(fid2,'%s\t%s\t%s\t%s\t%s\t%s\t\n',... 'DR ','DRds ','DRdm ','SNRnoi','SNR-ds','SNR-dm'); for j=1:nc

DBD



fprintf(fid2,'%4.2f\t%4.2f\t%4.2f\t%4.2f\t%4.2f\t%4.2f\t\n',... dr(j),drds(j),drdm(j),snr(j),snrds(j),snrdm(j)); end fprintf(fid2,'%4.2f\t%4.2f\t%4.2f\t%4.2f\t%4.2f\t%4.2f\t\n',... mdr,mdrds,mdrdm,msnr,msnrds,msnrdm); fclose(fid2); cd ..\..\..\algs2; %Plota e grava perfil medio, para posterior analise geral mpdp=20*log10(mean(10.^(pdp/20))); mpdpds=20*log10(mean(10.^(pdpds/20))); mpdpdm=20*log10(mean(10.^(pdpdm/20))); minimom=min(mpdp); if (minimom>-10) eixominm=-10; else eixominm=10*floor(minimom/10); end figure; plot (retardo,mpdp,'c');hold on; plot(retardo,mpdpds,':b','LineWidth',2); plot(retardo,mpdpdm,'r','LineWidth',2); title(sprintf('PDP medio para rota %s',arq)); xlabel('\tau [us]'); ylabel('P [dB]');grid on; axis([0 50 eixominm 0]); legend('STR','SLN','MLN');pause; arqg=sprintf('%s%smed',arq,parden); %grava fig media em jpg cd ..\results2\MF\figs; print('-dbmp',arqg); cd ..\..\..\algs2; hold off; disp('Salvando dados para analise .. ..'); destino=sprintf('..%sresults2%sMF%s%s%s','\','\','\',arq,parden); %nome do arqv para analise estatistica f1='pdp'; f2='pdpds'; f3='pdpdm'; save(destino, f1, f2, f3); %salva o resultado else disp('rota imprestavel!'); end

%Captura dados tese do Eduardo - SONDA STDCC clear; clc; cwd = pwd;cd(tempdir);pack;cd(cwd); %limpa tudo addpath ..\origem\STDCC\CETUC\andando; addpath ..\origem\STDCC\CETUC\fixo\; addpath ..\origem\STDCC\Copacabana\; disp ('Programa de Captura dos dados da tese do Eduardo'); arq = input('Nome do Arquivo : ','s'); disp('Carregando arquivo...'); arqn= sprintf('%s.txt',arq); fid=fopen(arqn); N=4096; %nr de pontos por varredura na tecnica por filtro casado yref=fscanf(fid,'%u',[1 ]); %yreference yinc=fscanf(fid,'%lg',[1 1]); %yincrement yo=fscanf(fid,'%lg',[1 1]); %yorigin c=1; try while (ftell(fid)~=-1) disp (sprintf('Varredura %u',c)); T(c)=fscanf(fid,'%lg',[1 1]); for i=1:N D(c,i)=fscanf(fid,'%u',[1 1]); Virg=fscanf(fid,'%1s',[1 1]); end c=c+1; clc; end catch disp('Fim do arquivo'); c=c-1; end fclose(fid); %--------------------------------------- %Trunca em 2555 pontos - Vide Tese do Edu - Ap A.3 pg 181-182 N2=2555; D=D(:,1:N2); %-------------------------------------------------------------------------------------- %Coverte valores em tensao disp('Convertendo valores lidos em tensao...'); pdpV=yo+(D-yref)*yinc; clear D; %------------------------------------- %Le manualmente Reference Level - precisa consultar arquivo *.rel RL(1:c)=0; disp(sprintf('Nr de varreduras: %d',c)); limitei=1;limitef=1; contador=1;

DBD



while limitef~=c RLin=input(sprintf('RL %d: ',contador)); limitef=input(sprintf('nr da varredura ateh onde RL %d vale: ',contador)); RL(limitei:limitef)=RLin; limitei=limitef+1; contador=contador+1; end disp('Reference Level atualizado'); %------------------------------------- %Converte valores de tensao em dB - Vide Tese do Edu - Ap A.3 pg 183 - Calibracao da tela pdpdb=78.82*pdpV-80; for i=1:c pdpdb(i,:)=pdpdb(i,:)-RL(i); end pdp=10.^(pdpdb/20); clear pdpdb pdpV; % Define vetor de retardo: resposta com 51,1 us em N2 ptos retardo=((1:N2)-1)*(51.1/N2); disp('PDP original em esc linear atualizado'); %------------------------------------- % Define parametros do denoising sh='s'; onda='sym8';COMP=''; esqlimiar='sqtwolog'; L=6; %OU L=7 parden=sprintf('%sL%d%s',sh,L,COMP); %------------------------------------- cwd = pwd;cd(tempdir);pack;cd(cwd); %limpa tudo disp('Realizando denoising...'); figure; jc=1;nc=0; for j=1:c pdpds(j,:) = abs(wden(pdp(j,:),esqlimiar,sh,'sln',L,onda)); pdpdm(j,:) = abs(wden(pdp(j,:),esqlimiar,sh,'mln',L,onda)); minimo=20*log10(min(pdp(j,:))); if (minimo>-10) eixomin=-10; else eixomin=10*floor(minimo/10); end maximo=20*log10(max(pdp(j,:))); eixomax=10*ceil(maximo/10); %Poe em escala log e plota pdpds(j,:)=20*log10(pdpds(j,:)); pdpdm(j,:)=20*log10(pdpdm(j,:)); pdp(j,:)=20*log10(pdp(j,:)); plot (retardo,pdp(j,:),'c');hold on; plot(retardo,pdpdm(j,:),'r','LineWidth',2); plot(retardo,pdpds(j,:),':b','LineWidth',2); title(sprintf('PDP %u de %u -- rota %s',j,c,arq)); xlabel('\tau [us]'); ylabel('P [dB]');grid on; legend('STR','MLN','SLN'); axis([0 50 eixomin eixomax]); hold off; ask=input('valida? ','s'); if (ask~='n')&(ask~='N') pdpc(jc,:)=pdp(j,:);pdpdsc(jc,:)=pdpds(j,:);pdpdmc(jc,:)=pdpdm(j,:); jc=jc+1;nc=jc-1; end end disp([int2str(c-nc) ' respostas cortadas']);

DBD


Apêndice B Algoritmos utilizados na estimação de AOA

Como no apêndice anterior, dada a semelhança dos algoritmos utilizados,

apenas alguns são apresentados aqui. O primeiro deles corresponde ao primeiro

bloco mencionado no capítulo 7, e é o responsável por gerar e armazenar os

vetores 3D de respostas temporais. % gera arquivo do tipo AOA--.mat contendo M (nr de sensores), % resolucao, e pdp_comp3d(t,rt,m) - o vetor 3D a ser processado para a analise de AOA nos algoritmos do tipo "analise" % 08/01/2003 % ------------------------------------------------------------------------------------- %Carga dos arquivos clear; clc; %limpa tudo addpath ..\funcoes\; addpath ..\resultados\; addpath ..\origem\leme; addpath ..\origem\gavea; disp ('Programa de Processamento dos Dados Medidos - direto'); disp ('Medidas paradas'); disp ('Sem tratamento de ruido'); disp (' '); disp (' Transforma do dominio da freqüência para o dominio dos retardos'); arq = input('Nome da Medida: ','s'); M= input('Nr de Sensores: '); span= 200; for m=1:M disp(sprintf('Carregando arquivos do sensor %g de %g ...', m,M)); arqm= sprintf('%s%dM.txt',arq,m); magn=load(arqm); arqf=sprintf('%s%dF.txt',arq,m); fase=load(arqf); clear arqm arqf; arqt=sprintf('%s%dT.txt',arq,m); tempo=load(arqt); clear arqm arqf arqt; disp('Arquivos carregados.'); [varreduras, pontos]=size(magn); %detecta tamanho pontos=800; %compatibiliza com tamanho do denoising % ------------------------CALCULO DOS PDPS---------------------------------------------------- magn=10.^(magn/20); for j=1:varreduras %compoe o numero complexo for p=1:pontos z(j,p)=magn(j,p)*cos(fase(j,p)*pi/180)+i*magn(j,p)*sin(fase(j,p)*pi/180); end end clear magn fase; w=BHarrisW(pontos); for j=1:varreduras %aplica a funcao janela zj(j,:)=z(j,:).*w; end clear z w; for j=1:varreduras aux2=ccorrelogramse(zj(j,:),pontos); aux2=fliplr(aux2'); aux3=angle(aux2); fase(j,:)=aux3(1,:); result(j,:)=abs(aux2); %Armazena como campo em vez de potencia, pois a estimacao espacial ira faze-lo clear aux2 aux3 ; end resolucao=10^3/span; limiret=1000; RT=round(limiret/resolucao); %--------------------------VETOR 3D PDP_COMP3D--------------------------------------------- for t=1:varreduras %compoe o numero complexo for rt=1:RT pdp_comp3d(t,rt,m)=result(t,rt)*cos(fase(t,rt))+i*result(t,rt)*sin(fase(t,rt)); end end vetor_T(m)=varreduras; % os tempos de varreduras sao desiguais para cada sensor

DBD


Apêndice B Algoritmos utilizados na estimação de AOA 282

clear fase result; end T=min(vetor_T); pdp_comp3d=pdp_comp3d(1:T,:,:); % ------SALVANDO DADOS PARA ANALISE ESTATISTICA - opcao AOA-------------------- disp('Salvando dados para analise...'); destino=sprintf('..%sresultados%s%sstr_stm','\','\',arq); f1='pdp_comp3d'; f2='resolucao'; save(destino, f1, f2); %salva o resultado clear destino; disp('Dados salvos.'); disp('Concluido.');clear;

O próximo algoritmo realiza a carga do vetor 3D especificado como

parâmetro de entrada, e estima o espectro espacial usando conformação de feixe. % Estimacao de espectro angular de arranjo ULA por beamforming % espacamento de 5 cm - fcentral = 1800 M - span 200 M %------------------------------------------------------------------ %Carga do arquivo clear; clc; %limpa tudo addpath ..\funcoes\; addpath ..\resultados\; addpath ..\origem\leme; addpath ..\origem\gavea; disp ('AOA - ULA - Abert Sint - beam');disp (' '); arq1=input('Nome do arquivo: ','s'); arq = input('Nome da extensao (str ou den):','s'); if arq=='den' arq='dpsmlnsym8L5'; end disp('Carregando arquivos...'); arq2= sprintf('%s%s_stm.mat',arq1,arq); load(arq2); %dados carregadas: pdp_comp3d(t,rt,m) e resolucao [T, RT, M]=size(pdp_comp3d); %detecta tamanho - T=tempo, RT=retardo, M=nr de sensores clear arq2; disp('Arquivo carregado...'); L=M*5; %nr de amostras no dominio AOA a serem observadas [-pi/2,pi/2] for a=1:L ang(a)=-pi/2+(a-1)*pi/(L-1); end angg=ang*180/pi; %vetor de AOA em graus ret=(1:RT)*resolucao; %define vetor de retardos [ns] for rt=1:RT for t=1:T for m=1:M Y(m,t)=pdp_comp3d(t,rt,m); end end %Y=flipud(Y); %Testa a direcao contraria de aquisicao no dominio espacial AOArt=beamform(Y,L,0.3); clear Y; %0.3 foi a rel de separacao entre dois elementos do array/compr onda (max = 0.5) for a=1:L AOA(rt,a)=AOArt(a); end clear AOArt; end disp('Espectros estimados.. ..'); AOA=10*log10(AOA); %converte valores de intensidade em escala log maxc=max(AOA,[],2);[maxa rtmax]=max(maxc); AOA=AOA-maxa; %grafico normalizado ao maximo npiso=min(min(AOA)); %novo piso figure;grid; mesh(angg,ret,AOA); xlabel('AOA [o]');ylabel('\tau [ns]'); zlabel('P [dB]'); title(sprintf('AOA x PDP: %s%s - Beamforming',arq1,arq)); arq3=sprintf('%s%sbeam3d.fig',arq1,arq); arq3b=sprintf('%s%sbeam3d',arq1,arq); arq3p=sprintf('%s%sbeam3dp.fig',arq1,arq); arq3bp=sprintf('%s%sbeam3dp',arq1,arq); cd ..\resultados\figuras\beam; saveas(gcf,arq3); print('-dbmp',arq3b); [th,r]=meshgrid(ang,ret); [retx,rety]=pol2cart(th,r); figure; mesh(retx,rety,AOA); title(sprintf('AOA x PDP: %s%s - Beamforming',arq1,arq)); xlabel('\tau [ns]');ylabel('\tau [ns]'); zlabel('P [dB]'); text(ret(RT-5),-ret(RT),max(AOA(RT,:))+10,'-\pi/2'); text(ret(RT-5),ret(RT),max(AOA(RT,:))+10,'\pi/2'); saveas(gcf,arq3p); print('-dbmp',arq3bp); cd ..\..\..\algoritmos; %Analise do nr de picos e determinacao dos picos for rt=1:RT picos=0; for a=2:L-1 ant=AOA(rt,a-1); pos=AOA(rt,a+1); if (AOA(rt,a)>ant)&(AOA(rt,a)>pos) picos=picos+1; AOAd_a(rt,picos)=angg(a);

DBD



AOAd_p(rt,picos)=AOA(rt,a); %AOAd(rt,a)=AOA(rt,a); else %AOAd(rt,a)=npiso; end end %AOAd(rt,1)=0; %AOAd(rt,L)=0; picostot(rt)=picos; end %Mostra dominio AOA, retardo a retardo, destacando os picos figure; cd ..\resultados\figuras\beam; for rt=(rtmax-5):(rtmax+5) plot(angg,AOA(rt,:),'LineWidth',2); grid; axis([-90 90 -60 0]); xlabel('AOA [o]');ylabel('P [dB]'); title(['AOA x PDP (\tau = ',int2str(rt*resolucao),' ns) - ',arq1,arq,' - Beamforming']); hold on; for p=1:picostot(rt) %Busca picos do retardo rt para destaca-los no grafico plot(AOAd_a(rt,p),AOAd_p(rt,p),'d','LineWidth',2); end hold off; arq4=sprintf('%s%sbeamrt%d.fig',arq1,arq,rt*resolucao); saveas(gcf,arq4); arq4b=sprintf('%s%sbeamrt%d',arq1,arq,rt*resolucao); print('-djpeg',arq4b); end cd ..\..\..\algoritmos; %grava em arquivo texto os resultados arq5=sprintf('%s%sbeamrt%d.xls',arq1,arq,rtmax*resolucao); cd ..\resultados\planilhas\beam; fid = fopen(arq5,'w'); for rt=(rtmax-5):(rtmax+5) fprintf(fid,'\n%s\t','rt'); fprintf(fid,'%d\n',rt*resolucao); fprintf(fid,'%s\t','AOA [o]'); for p=1:picostot(rt) fprintf(fid,'%4.2f\t',AOAd_a(rt,p)); end fprintf(fid,'\n%s\t','P [dB]'); for p=1:picostot(rt) fprintf(fid,'%4.2f\t',AOAd_p(rt,p)); end end fclose(fid); cd ..\..\..\algoritmos;

Os algoritmos para os demais métodos espaciais utilizados são apresentados

a seguir, em versão reduzida, pois a apresentação dos gráficos com os espectros e

o armazenamento de dados para análise posterior são análogos. % Estimacao de espectro angular de arranjo ULA por Capon %------------------------------------------------------------------ %Carga do arquivo clear; clc; %limpa tudo addpath ..\funcoes\; addpath ..\resultados\; addpath ..\origem\leme; addpath ..\origem\gavea; disp ('AOA - ULA - Abert Sint - CAPON'); disp (' '); arq1=input('Nome do arquivo: ','s'); arq = input('Nome da extensao (str ou den):','s'); if arq=='den' arq='dpsmlnsym8L5'; end disp('Carregando arquivos...'); arq2= sprintf('%s%s_stm.mat',arq1,arq); load(arq2); %dados carregadas: pdp_comp3d(t,rt,m) e resolucao [T, RT, M]=size(pdp_comp3d); %detecta tamanho - T=tempo, RT=retardo, M=nr de sensores clear arq2; disp('Arquivo carregado...'); L=M*5; %nr de amostras no dominio AOA a serem observadas [-pi/2,pi/2] for a=1:L ang(a)=-pi/2+(a-1)*pi/(L-1); end angg=ang*180/pi; %vetor de AOA em graus ret=(1:RT)*resolucao; %define vetor de retardos [ns] for rt=1:RT for t=1:T for m=1:M Y(m,t)=pdp_comp3d(t,rt,m);

DBD



end end AOArt=capon_sp(Y,L,0.3); clear Y; for a=1:L AOA(rt,a)=AOArt(a); end clear AOArt; end disp('Espectros estimados.. ..'); AOA=10*log10(AOA); %converte valores de intensidade em escala log

% Estimacao de espectro angular de arranjo ULA por Spectral MUSIC, com Nr ftes estimado pelo nr de picos do beamf. %------------------------------------------------------------------ %Carga do arquivo clear; clc; %limpa tudo addpath ..\funcoes\; addpath ..\resultados\; addpath ..\origem\leme; addpath ..\origem\gavea; disp ('AOA - ULA - Abert Sint - beam + MUSIC'); disp (' '); arq1=input('Nome do arquivo: ','s'); arq = input('Nome da extensao (str ou den):','s'); if arq=='den' arq='dpsmlnsym8L5'; end disp('Carregando arquivos...'); arq2= sprintf('%s%s_stm.mat',arq1,arq); load(arq2); %dados carregadas: pdp_comp3d(t,rt,m) e resolucao [T, RT, M]=size(pdp_comp3d); %detecta tamanho - T=tempo, RT=retardo, M=nr de sensores clear arq2; disp('Arquivo carregado...'); L=M*5; %nr de amostras no dominio AOA a serem observadas [-pi/2,pi/2] for a=1:L ang(a)=-pi/2+(a-1)*pi/(L-1); end angg=ang*180/pi; %vetor de AOA em graus ret=(1:RT)*resolucao; %define vetor de retardos [ns] for rt=1:RT for t=1:T for m=1:M Y(m,t)=pdp_comp3d(t,rt,m); end end AOArt=beamform(Y,L,0.3); clear Y; for a=1:L AOA(rt,a)=AOArt(a); end clear AOArt; end disp('Espectros por conf. de feixe estimados.. ..'); AOA=10*log10(AOA); %converte valores de intensidade em escala log maxc=max(AOA,[],2);[maxp rtmax]=max(maxc); AOA=AOA-maxp; %grafico normalizado ao maximo npiso=min(min(AOA)); %novo piso %Analise do nr de picos e determinacao dos picos for rt=1:RT picos=0; for a=2:L-1 ant=AOA(rt,a-1); pos=AOA(rt,a+1); if (AOA(rt,a)>ant)&(AOA(rt,a)>pos) picos=picos+1; end end picostot(rt)=picos; %Estima espectro MUSIC, considerando como nr de picos o vetor estimado a partir da estimativa por beamf., picostot for t=1:T for m=1:M Y(m,t)=pdp_comp3d(t,rt,m); end end np=picostot(rt); ampert=mean(mean(abs(Y))); %cria valor ampert, com amplit estim. do pico de retardo rt para pond. o pseudo-espec. AOArt=spectral_music_doa(Y,L,np,0.3); clear Y; for a=1:L AOAmusic(rt,a)=AOArt(a)*ampert^2; end picosm=0; for a=2:L-1 ant=AOAmusic(rt,a-1); pos=AOAmusic(rt,a+1);

DBD



if (AOAmusic(rt,a)>ant)&(AOAmusic(rt,a)>pos) picosm=picosm+1; AOAmd_a(rt,picosm)=angg(a); AOAmd_p(rt,picosm)=AOAmusic(rt,a); end end picostotm(rt)=picosm; clear AOArt ampert; end clear AOA; disp('Espectros MUSIC estimados!'); AOA=10*log10(abs(AOAmusic)); %converte valores de intensidade em escala log

% Estimacao de espectro angular de arranjo ULA por ESPRIT, com Nr ftes estimado pelo nr de picos do beamforming %------------------------------------------------------------------ %Carga do arquivo clear; clc; %limpa tudo addpath ..\funcoes\; addpath ..\resultados\; addpath ..\origem\leme; addpath ..\origem\gavea; disp ('AOA - ULA - Abert Sint - beam + ESPRIT'); disp (' '); arq1=input('Nome do arquivo: ','s'); arq = input('Nome da extensao (str ou den):','s'); if arq=='den' arq='dpsmlnsym8L5'; end disp('Carregando arquivos...'); arq2= sprintf('%s%s_stm.mat',arq1,arq); load(arq2); %dados carregadas: pdp_comp3d(t,rt,m) e resolucao [T, RT, M]=size(pdp_comp3d); %detecta tamanho - T=tempo, RT=retardo, M=nr de sensores clear arq2; disp('Arquivo carregado...'); L=M*5; %nr de amostras no dominio AOA a serem observadas [-pi/2,pi/2] for a=1:L ang(a)=-pi/2+(a-1)*pi/(L-1); end angg=ang*180/pi; %vetor de AOA em graus ret=(1:RT)*resolucao; %define vetor de retardos [ns] for rt=1:RT for t=1:T for m=1:M Y(m,t)=pdp_comp3d(t,rt,m); end end AOArt=beamform(Y,L,0.3); clear Y; for a=1:L AOA(rt,a)=AOArt(a); end clear AOArt; end disp('Espectros por conf. de feixe estimados.. ..'); AOA=10*log10(AOA); %converte valores de intensidade em escala log maxc=max(AOA,[],2);[maxp rtmax]=max(maxc); AOA=AOA-maxp; %grafico normalizado ao maximo npiso=min(min(AOA)); %novo piso %Analise do nr de picos e determinacao dos picos AOAesprit(1:RT,1:L)=0; for rt=1:RT picos=0; for a=2:L-1 ant=AOA(rt,a-1); pos=AOA(rt,a+1); if (AOA(rt,a)>ant)&(AOA(rt,a)>pos) picos=picos+1; end end picostot(rt)=picos; %Estima DOAs por ESPRIT, considerando como nr de picos o vetor estimado a partir da estimativa por beamf, picostot for t=1:T for m=1:M Y(m,t)=pdp_comp3d(t,rt,m); end end ampe(rt)=mean(mean(abs(Y))); %cria ampert, com amplit estim. do pico de retardo rt para ponderar o pseudo-espectro np=picostot(rt); AOAe=esprit_doa(Y,np,0.3)'; clear Y; AOAe=real(AOAe);

DBD



% A rotina abaixo mapeia apenas as DOAs validas (~= 90) para a matriz final AOAesprit, % que tem L-1 bins angulares, e atualiza o vetor picostote! picostote(rt)=0; for a=1:np if abs(AOAe(a))~=90 i=round((L-1)*(AOAe(a)+90)/180)+1; AOAesprit(rt,i)=ampe(rt)^2; picostote(rt)=picostote(rt)+1; AOAd_a(rt,picostote(rt))=AOAe(a); end end end clear AOA; %Normaliza pelo maximo maxa=max(ampe)^2; mina=min(ampe)^2; AOAesprit=AOAesprit/maxa; mina=mina/maxa; %Mapeia os valores nulos da matriz AOAesprit, que correspondem na verdade a valores de inicializaçao, %em valores iguais ao piso (mina) for rt=1:RT for a=1:L if (AOAesprit(rt,a)<mina) AOAesprit(rt,a)=mina; end end end disp('Espectros ESPRIT estimados!'); AOA=10*log10(AOAesprit); %converte valores de intensidade em escala log

Os algoritmos de estimação do espectro espacial utilizam funções acessadas

externamente, cujos códigos vêm a seguir. Como no apêndice anterior, as funções

abaixo foram extraídas de [16], mas incorporam algumas modificações que as

tornam mais adequadas aos algoritmos principais utilizados. function phi=beamform(Y,L,d) % phi=beamform(Y,L,d); % Y <- dados do arranjo ULA % L <- nr de amostras em [-pi/2,pi/2] para procurar pelas fontes % d <- espaçamento entre sensores em comprimentos de onda % phi <- estimativas espectrais em L ângulos igualmente espaçados (L-1 na verdade - em 13/01/2003) em [-90,90] graus [m,N]=size(Y); % computa a matriz covariância amostral R=Y*Y'/N; phi=zeros(L,1); for i = 1 : L, a=exp(-2*pi*j*d*sin(-pi/2 + pi*(i-1)/(L-1))*[0:m-1].'); phi(i)=real(a'*R*a); end %OBS: Na linha 24, modifiquei "pi*(i-1)/L" para "pi*(i-1)/(L-1)" em 13/01/2003

function phi=capon_sp(Y,L,d) % phi=capon_sp(Y,L,d); % Y <- dados do arranjo ULA % L <- nr de amostras em [-pi/2,pi/2] para procurar pelas fontes % d <- espaçamento entre sensores em comprimentos de onda % phi <- estimativas espectrais em L ângulos igualmente espaçados (L-1 na verdade - em 13/01/2003) em [-90,90] graus [m,N]=size(Y); % computa a inversa da matriz covariância amostral %IR=inv(Y*Y'/N); IR=(1/N)*(Y*Y')\eye(m); %Esta istrucao eh analoga e mais rapida que a anterior phi=zeros(L,1); for i = 1 : L, a=exp(-2*pi*j*d*sin(-pi/2 + pi*(i-1)/(L-1))*[0:m-1].'); phi(i)=1/real(a'*IR*a); end %OBS: Na linha 24, modifiquei "pi*(i-1)/L" para "pi*(i-1)/(L-1)" em 13/01/2003

function phi=spectral_music_doa(Y,L,n,d) % phi=spectral_music_doa(Y,L,n,d)

DBD



% Y <- dados do arranjo ULA % L <- nr de amostras em [-pi/2,pi/2] para procurar pelas fontes % n <- nr de fontes % d <- espaçamento entre sensores em comprimentos de onda % phi <- estimativas espectrais em L ângulos igualmente espaçados (L-1 na verdade - em 13/01/2003) em [-90,90] graus [m,N]=size(Y); % computa a matriz covariância amostral R=Y*Y'/N; % realiza a autodecomposição, usando SVD pois ele ordena os autovalores [U,D,V]=svd(R); G=U(:,n+1:m); GG = G*G'; phi=zeros(L,1); for i = 1 : L a=exp(-2*pi*j*d*sin(-pi/2 + pi*(i-1)/(L-1))*[0:m-1].'); phi(i)=real(inv(a'*GG*a)); end %OBS: Na linha 29, modifiquei "pi*(i-1)/L" para "pi*(i-1)/(L-1)" em 13/01/2003

function doa=esprit_doa(Y,n,d) % doa=esprit_doa(Y,n,d) % Y <- dados do arranjo ULA % n <- nr de fontes % d <- espaçamento entre sensores em comprimentos de onda % doa -> vetor com as estimativas de DOA [m,N]=size(Y); % computa a matriz covariância amostral R=Y*Y'/N; % realiza a autodecomposição, usando SVD pois ele ordena os autovalores [U,D,V]=svd(R); S=U(:,1:n); phi = S(1:m-1,:)\S(2:m,:); w=-angle(eig(phi)); doa=asin(w/d/pi/2)*180/pi;

Por fim, o algoritmo correspondente ao primeiro módulo do bloco de

processamento das medidas com o arranjo real é apresentado. As rotinas de

análise do capítulo 8 são pratiamente as mesmas do capiítulo 7. % arquivo lido do tipo "leme I J V 801 K.txt % onde I eh a posicao ao longo da rota (nr) % J eh a disp relativa ('n' ou 'z') % V eh a varredura (30 para cetuc1 e 3 para todas as demais) % K eh 'F', 'M' ou 'T', como na estrutura original da sonda do LH % Este aqui roda direto, gerando os arquivos necess'arios de uma so vez % gera arq do tipo nome--.mat contendo resolucao, e pdp_comp3d(t,rt,m) - o vetor 3D a ser processado na analise de AOA % 04/02/2003 % ------------------------------------------------------------------------------------- %Carga dos arquivos clear; clc; %limpa tudo addpath ..\funcoes\; addpath ..\resultados\; addpath ..\origem\leme; addpath ..\origem\cetuc; disp ('Programa de organizacao dos dados medidos'); disp ('ULA16'); disp ('Transforma do dominio da freqüência para o dominio dos retardos'); disp ('Le todos os pontos da rota'); arq = input('Nome da Medida (leme - L ou cetuc - C): ','s'); nz = input('N ou Z:','s'); if (arq=='L')|(arq=='l')|(arq=='leme') np= 10; arq2='leme'; elseif (arq=='C')|(arq=='c')|(arq=='cetuc') np= 5; arq2='cetuc'; else np=0; end M=16; %ULA 16 pontos=801; span= 200; resolucao=10^3/span; %nano segundos limiret=2000; RT=round(limiret/resolucao); w=BHarrisW(pontos); %funcao janela for n=1:np if (np==5)&(n==1) nv=30; else nv=3; end

DBD



disp(sprintf('Carregando arquivos da posicao %g de %g ...', n,np)); for v=1:nv arqm= sprintf('%s%d%s%d801M.txt',arq2,n,nz,v); magn=load(arqm); magn=magn.'; %transposta nao conjugada magn3d(v,:,:)=magn; clear magn; arqf=sprintf('%s%d%s%d801F.txt',arq2,n,nz,v); fase=load(arqf); fase=fase.'; fase3d(v,:,:)=fase; clear fase; arqt=sprintf('%s%d%s%d801T.txt',arq2,n,nz,v); tempo=load(arqt); tempo=tempo.'; tempo3d(v,:,:)=tempo; clear tempo; clear arqm arqf arqt; %estrutura dos arqs tempo3d,fase3d e magn3d: (v,rt,m) end disp('Arquivos carregados.'); clear magn fase tempo; % ------------------------CALCULO DOS PDPS---------------------------------------------------- for m=1:M for v=1:nv magn(v,:)=10.^(magn3d(v,:,m)/20); fase(v,:)=fase3d(v,:,m); zj(v,:)=(magn(v,:).*cos(fase(v,:)*pi/180)+i*magn(v,:).*sin(fase(v,:)*pi/180)).*w; aux2=ccorrelogramse(zj(v,:), pontos); %aplica o correlograma aux2=fliplr(aux2'); aux3=angle(aux2); fase(v,:)=aux3(1,:); %GUARDA RESPOSTA DE FASE result(v,:)=abs(aux2); %Armazena como campo em vez de potencia, pois a estimacao espacial ira faze-lo clear aux2 aux3 ; %--------------------------VETOR 3D PDP_COMP3D TRUNCADO------------------------------ for rt=1:RT pdp_comp3d(v,rt,m)=result(v,rt)*cos(fase(v,rt))+i*result(v,rt)*sin(fase(v,rt)); end end disp (sprintf('PDPs calculados para o sensor %d',m)); clear magn zj fase result; end disp(sprintf('Salvando dados para %s%d%s',arq2,n,nz)); destino=sprintf('..%sresultados%s%s%d%s_stm','\','\',arq2,n,nz); f1='pdp_comp3d'; f2='resolucao'; save(destino, f1, f2); %salva o resultado clear destino pdp_comp3d magn3d fase3d tempo3d; end disp('Dados salvos.'); disp('Concluido.');

DBD


Documents

10 Referências bibliográficas - PUC-Rio · 10 Referências bibliográficas 1 FUHL, J. Smart antennas for second and third generation mobile communications systems. Viena - Áustria,