Inverse of Matrix and Signal reconstruction
in [Matlab]
|
From: kk KKsingh on 22 Jul 2010 16:47 Hi ! I have a matrix which looks like (A*A) where A is Fourier Matrix,It need to apply in this form (A* WA)^(inverse) * A*Wy; Where A* is conjugate of A....y is irregular sample signal (So A*Wy; is bandlimited spectra)...W is the weights for samples... which is the best way to invert is using backslash (\) is making the condition number bad resulting in unstable result. same happens with using inv(A*W*A)...What you suggest about pinv(A*WA).... Also did any body tried writing this code for signal processing, I wasted so much time on this algorithm still not able to get any results ...Following is my algorithm 1. Take a irregular signal Y 2.Apply A*Y; (Make it band limited , cut off high frequencies) 3. Apply inv(A*WA) to A*Y; 4. Apply IDFT to get samples back Please let me know, if any one has ever tried this ..I really need some help on this
From: Walter Roberson on 22 Jul 2010 17:48 kk KKsingh wrote: > I have a matrix which looks like (A*A) where A is Fourier Matrix,It need > to apply in this form > > (A* WA)^(inverse) * A*Wy; > > Where A* is conjugate of A Is it really the conjugate of A, conj(A) ? Or is it the conjugate transpose of A, A' ? If it is conj(A) then unless A is square you have a dimension problem as the second dimension of conj(A) would not match the first dimension of mult(W,A) . .....y is irregular sample signal (So A*Wy; is > bandlimited spectra)...W is the weights for samples... > > which is the best way to invert is using backslash (\) is making the > condition number bad resulting in unstable result. same happens with > using inv(A*W*A)...What you suggest about pinv(A*WA).... The pseudo-inverse of a matrix is not unique, but pinv() choses one particular value of it and \ chooses a different value of it (see the pinv documentation.) If your data is so badly conditioned then it is not immediately clear to me that either of the two matrices would be the "right" matrices for your purposes (but I would need more research to be certain that one or the other would not be meaningful for your problem.) What I _suspect_ this bad condition number is telling you is that your model of the process is inadequate. > Also did any body tried writing this code for signal processing, I > wasted so much time on this algorithm still not able to get any results > ...Following is my algorithm > > 1. Take a irregular signal Y > 2.Apply A*Y; (Make it band limited , cut off high frequencies) > 3. Apply inv(A*WA) to A*Y; > 4. Apply IDFT to get samples back > Please let me know, if any one has ever tried this .. Yes, that kind of approach has been tried (though a transpose might perhaps be needed in there somewhere.) I do not know the conditions under which it can be successful.
From: kk KKsingh on 22 Jul 2010 18:10 Walter Roberson <roberson(a)hushmail.com> wrote in message <i2aein$rj$1(a)canopus.cc.umanitoba.ca>... > kk KKsingh wrote: > > > I have a matrix which looks like (A*A) where A is Fourier Matrix,It need > > to apply in this form > > > > (A* WA)^(inverse) * A*Wy; > > > > Where A* is conjugate of A > > Is it really the conjugate of A, conj(A) ? Or is it the conjugate transpose > of A, A' ? If it is conj(A) then unless A is square you have a dimension > problem as the second dimension of conj(A) would not match the first dimension > of mult(W,A) . > > > ....y is irregular sample signal (So A*Wy; is > > bandlimited spectra)...W is the weights for samples... > > > > which is the best way to invert is using backslash (\) is making the > > condition number bad resulting in unstable result. same happens with > > using inv(A*W*A)...What you suggest about pinv(A*WA).... > > The pseudo-inverse of a matrix is not unique, but pinv() choses one particular > value of it and \ chooses a different value of it (see the pinv > documentation.) If your data is so badly conditioned then it is not > immediately clear to me that either of the two matrices would be the "right" > matrices for your purposes (but I would need more research to be certain that > one or the other would not be meaningful for your problem.) > > What I _suspect_ this bad condition number is telling you is that your model > of the process is inadequate. > > > > Also did any body tried writing this code for signal processing, I > > wasted so much time on this algorithm still not able to get any results > > ...Following is my algorithm > > > > 1. Take a irregular signal Y > > 2.Apply A*Y; (Make it band limited , cut off high frequencies) > > 3. Apply inv(A*WA) to A*Y; > > 4. Apply IDFT to get samples back > > > Please let me know, if any one has ever tried this .. > > Yes, that kind of approach has been tried (though a transpose might perhaps be > needed in there somewhere.) I do not know the conditions under which it can be > successful. Its Conjugate Transpose
From: Walter Roberson on 22 Jul 2010 18:48 kk KKsingh wrote: > Walter Roberson <roberson(a)hushmail.com> wrote in message > <i2aein$rj$1(a)canopus.cc.umanitoba.ca>... >> kk KKsingh wrote: >> >> > I have a matrix which looks like (A*A) where A is Fourier Matrix,It >> need > to apply in this form >> > > (A* WA)^(inverse) * A*Wy; >> > > Where A* is conjugate of A > Its Conjugate Transpose Is the second dimension of A greater than the first dimension? If so then (A'*A) will not have an inverse, it appears. You mention that A is a Fourier Matrix, but other parts of your formulae become more natural to understand if A is a row vector rather than a matrix. Is A a 2D matrix? If so, then how were the fourier coefficients determined? For it to be a 2D fourier matrix would be suggestive of a 2D Fourier transform, but your mention of IDFT suggests that you are working in 1D not 2D.
From: kk KKsingh on 23 Jul 2010 00:49
Walter Roberson <roberson(a)hushmail.com> wrote in message <i2ai3k$61n$1(a)canopus.cc.umanitoba.ca>... > kk KKsingh wrote: > > Walter Roberson <roberson(a)hushmail.com> wrote in message > > <i2aein$rj$1(a)canopus.cc.umanitoba.ca>... > >> kk KKsingh wrote: > >> > >> > I have a matrix which looks like (A*A) where A is Fourier Matrix,It > >> need > to apply in this form > >> > > (A* WA)^(inverse) * A*Wy; > >> > > Where A* is conjugate of A > > > Its Conjugate Transpose > > Is the second dimension of A greater than the first dimension? If so then > (A'*A) will not have an inverse, it appears. > > You mention that A is a Fourier Matrix, but other parts of your formulae > become more natural to understand if A is a row vector rather than a matrix. > Is A a 2D matrix? If so, then how were the fourier coefficients determined? > For it to be a 2D fourier matrix would be suggestive of a 2D Fourier > transform, but your mention of IDFT suggests that you are working in 1D not 2D. %This is to test codes clear all; close all; N=1024; %Number of Sample points Fs=1000;%Sampling Frequency dt=1/Fs;%Sampling interval fo=50;% Frequency of sin wave t=0:dt:dt*(N-1);%Sampling Period omega=2*pi*fo; S=sin(omega*t); %Applying FFT f=fftshift(fft(S)); %Applying the decimation kernal pd=input('Enter the percent decimation') if pd == 0 W = ones(size(S)); elseif pd ==1 W = zeros(size(S)); else W = rand(size(S)); W = W + (1-mean(W) - pd); mean(W) W = round(W); end S1=W.*S; %finding zeros from the signal index=find(S1==0) %Removing zeros from the signal S2=S1; t2=t; S2(index)=[]; t2(index)=[]; N1=numel(S2); freq_NDFT=Fs/N1*(-N/2:N/2-1); S2=S2(:); t2=t2(:); freq_NDFT=freq_NDFT(:); %Forward Fourier Transform DFT=exp(-2*pi*i*freq_NDFT*t2'); X_DFT=DFT*S2; %------------------------------------------------------ %Schonvillie NDFT is solved using Reinmann sum weights are %used here %NDFT Accoring to Schonvillie ยจ N3=numel(S2); Period=max(t2)-min(t2); deltax_n=zeros(N1,1); deltax_n(1,:)=((t2(2)-(t2(N1)-Period)))/2; deltax_n(N1,:)=((t2(1)+Period)-(t2(N3-1)))/2; for k=2:N1-1; deltax_n(k,:)= ((t2(k+1)-t2(k-1))/2); end %Now applying weights to the DFT X_DFTw=DFT*(S2.*deltax_n) %Making inverse kernal X_DFT_least= pinv(DFT*diag(deltax_n)*DFT')*X_DFTw; figure(1) subplot(311) plot(freq_NDFT,abs(X_DFT)) legend('NDFT') subplot(312) plot(freq_NDFT,abs(X_DFTw)) legend('NDFT with Weights') subplot(313) plot(freq_NDFT,abs(X_DFT_least)) legend('Least Squares NDFT') Please let me know if there is any mistake ! Except bandlimitation i have done every thing ! I will band limit it by setting high frequency to zero.....will be taking only 20% samples kk |