Hi All, I am working on a simulation in Matlab. Basically I have 12 random spread spectrum signals sampled. Then I have Matlab randomly pick one of these and add AWGN on. Finally, I use matched filtering of this random signal to each of the 12 predefined signals and repeat this several times at various SNRs. The goal is obviously to detect which of the 12 predefined signals is being used here. Basically for my matched filter in Matlab, I simply take the following: max(conv(x_r,fliplr(x1))) and repeat this for x_r (random signal) against all of the originals. The one with the maximum output should be the original signal that the random one was generated from. Unfortunately I am getting extremely high error results, between 60-70% with not real improvement even when SNR is increased significantly. I have verified the signals are looking very similar and have even had the same result when I take the noise out so obviously there is a major flaw with my filter. It is also worth mentioning that I do make sure to normalize all signals as well. If anyone has any suggestions of what I am doing wrong or another way to get better results soon it would be greatly appreciated. Thanks!
Matched Filtering using Matlab
Started by ●July 24, 2008
Reply by ●July 24, 20082008-07-24
On Jul 24, 11:25 am, "squick40" <steve.qu...@gmail.com> wrote:> Hi All, > > I am working on a simulation in Matlab. Basically I have 12 random spread > spectrum signals sampled. Then I have Matlab randomly pick one of these > and add AWGN on. Finally, I use matched filtering of this random signal to > each of the 12 predefined signals and repeat this several times at various > SNRs. The goal is obviously to detect which of the 12 predefined signals > is being used here. > > Basically for my matched filter in Matlab, I simply take the following: > > max(conv(x_r,fliplr(x1))) > > and repeat this for x_r (random signal) against all of the originals. The > one with the maximum output should be the original signal that the random > one was generated from. Unfortunately I am getting extremely high error > results, between 60-70% with not real improvement even when SNR is > increased significantly. I have verified the signals are looking very > similar and have even had the same result when I take the noise out so > obviously there is a major flaw with my filter. It is also worth > mentioning that I do make sure to normalize all signals as well. > > If anyone has any suggestions of what I am doing wrong or another way to > get better results soon it would be greatly appreciated. Thanks!Your approach sounds reasonable. How are you normalizing the signals? What do their cross-correlations look like? Ideally, you'd like all 12 signals to be perfectly orthogonal to each other such that you get zero at the correct sampling instant when you correlate any pair of two unlike signals. At the same time, you want the peak of the autocorrelation of each signal to be large compared to its sidelobes so you can accurately detect chip timing. Since what you're really doing is a cross-correlation between your "received signal" and your 12 template signals, you could replace the conv() with a call to xcorr() to take one piece of complexity (the time-flipping) out of the equation. In the noiseless case, as long as your signals are scaled properly and you assume a priori knowledge of the chip timing, you should not see any errors at all. Perhaps you could post your code. Jason
Reply by ●July 24, 20082008-07-24
>On Jul 24, 11:25 am, "squick40" <steve.qu...@gmail.com> wrote: >> Hi All, >> >> I am working on a simulation in Matlab. Basically I have 12 randomspread>> spectrum signals sampled. Then I have Matlab randomly pick one ofthese>> and add AWGN on. Finally, I use matched filtering of this randomsignal to>> each of the 12 predefined signals and repeat this several times atvarious>> SNRs. The goal is obviously to detect which of the 12 predefinedsignals>> is being used here. >> >> Basically for my matched filter in Matlab, I simply take thefollowing:>> >> max(conv(x_r,fliplr(x1))) >> >> and repeat this for x_r (random signal) against all of the originals.The>> one with the maximum output should be the original signal that therandom>> one was generated from. Unfortunately I am getting extremely higherror>> results, between 60-70% with not real improvement even when SNR is >> increased significantly. I have verified the signals are looking very >> similar and have even had the same result when I take the noise out so >> obviously there is a major flaw with my filter. It is also worth >> mentioning that I do make sure to normalize all signals as well. >> >> If anyone has any suggestions of what I am doing wrong or another wayto>> get better results soon it would be greatly appreciated. Thanks! > >Your approach sounds reasonable. How are you normalizing the signals? >What do their cross-correlations look like? Ideally, you'd like all 12 >signals to be perfectly orthogonal to each other such that you get >zero at the correct sampling instant when you correlate any pair of >two unlike signals. At the same time, you want the peak of the >autocorrelation of each signal to be large compared to its sidelobes >so you can accurately detect chip timing. Since what you're really >doing is a cross-correlation between your "received signal" and your >12 template signals, you could replace the conv() with a call to >xcorr() to take one piece of complexity (the time-flipping) out of the >equation. In the noiseless case, as long as your signals are scaled >properly and you assume a priori knowledge of the chip timing, you >should not see any errors at all. Perhaps you could post your code. > >Jason >Thanks Jason, Below is a copy of my code which is rather long. I have tried the xcorr() function with the same results. Let me know if you have any suggestions. I really do appreciate the help. snr = -20:10:20; % SNR values to be examined ns = length(snr) % Total number of SNR values to be used N = 128 % Spreading code length; T = NT_c; T_c = 1/W =chip duration Ns = 4 % Upsampling factor qo = sign(randn(N,1)); % Random set q = kron(qo,ones(Ns,1)); % Upsampled spreading code q = q/norm(q); % Normalize spreading code to unit energy % Channel Realizations L = 20 % normalized channel multipath spread at chip rate (L+1 multipath components) % Channel 11 Np11 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec11 = round(L*Ns*rand(Np11,1)); % The location of path delays at the upsampling res. tau_vec11 = tau_vec11 - min(tau_vec11)+1; % Normalized path delays h11 = zeros(L*Ns+1,1); % initialization of channel h11(tau_vec11) = sqrt(.5)*(randn(Np11,1) + j*randn(Np11,1)); % complex channel path gains % Channel 12 Np12 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec12 = round(L*Ns*rand(Np12,1)); % The location of path delays at the upsampling res. tau_vec12 = tau_vec12 - min(tau_vec12)+1; % Normalized path delays h12 = zeros(L*Ns+1,1); % initialization of channel h12(tau_vec12) = sqrt(.5)*(randn(Np12,1) + j*randn(Np12,1)); % complex channel path gains % Channel 13 Np13 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec13 = round(L*Ns*rand(Np13,1)); % The location of path delays at the upsampling res. tau_vec13 = tau_vec13 - min(tau_vec13)+1; % Normalized path delays h13 = zeros(L*Ns+1,1); % initialization of channel h13(tau_vec13) = sqrt(.5)*(randn(Np13,1) + j*randn(Np13,1)); % complex channel path gains % Channel 14 Np14 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec14 = round(L*Ns*rand(Np14,1)); % The location of path delays at the upsampling res. tau_vec14 = tau_vec14 - min(tau_vec14)+1; % Normalized path delays h14 = zeros(L*Ns+1,1); % initialization of channel h14(tau_vec14) = sqrt(.5)*(randn(Np14,1) + j*randn(Np14,1)); % complex channel path gains % Channel 21 Np21 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec21 = round(L*Ns*rand(Np21,1)); % The location of path delays at the upsampling res. tau_vec21 = tau_vec21 - min(tau_vec21)+1; % Normalized path delays h21 = zeros(L*Ns+1,1); % initialization of channel h21(tau_vec21) = sqrt(.5)*(randn(Np21,1) + j*randn(Np21,1)); % complex channel path gains % Channel 22 Np22 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec22 = round(L*Ns*rand(Np22,1)); % The location of path delays at the upsampling res. tau_vec22 = tau_vec22 - min(tau_vec22)+1; % Normalized path delays h22 = zeros(L*Ns+1,1); % initialization of channel h22(tau_vec22) = sqrt(.5)*(randn(Np22,1) + j*randn(Np22,1)); % complex channel path gains % Channel 23 Np23 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec23 = round(L*Ns*rand(Np23,1)); % The location of path delays at the upsampling res. tau_vec23 = tau_vec23 - min(tau_vec23)+1; % Normalized path delays h23 = zeros(L*Ns+1,1); % initialization of channel h23(tau_vec23) = sqrt(.5)*(randn(Np23,1) + j*randn(Np23,1)); % complex channel path gains % Channel 24 Np24 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec24 = round(L*Ns*rand(Np24,1)); % The location of path delays at the upsampling res. tau_vec24 = tau_vec24 - min(tau_vec24)+1; % Normalized path delays h24 = zeros(L*Ns+1,1); % initialization of channel h24(tau_vec24) = sqrt(.5)*(randn(Np24,1) + j*randn(Np24,1)); % complex channel path gains % Channel 31 Np31 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec31 = round(L*Ns*rand(Np31,1)); % The location of path delays at the upsampling res. tau_vec31 = tau_vec31 - min(tau_vec31)+1; % Normalized path delays h31 = zeros(L*Ns+1,1); % initialization of channel h31(tau_vec31) = sqrt(.5)*(randn(Np31,1) + j*randn(Np31,1)); % complex channel path gains % Channel 32 Np32 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec32 = round(L*Ns*rand(Np32,1)); % The location of path delays at the upsampling res. tau_vec32 = tau_vec32 - min(tau_vec32)+1; % Normalized path delays h32 = zeros(L*Ns+1,1); % initialization of channel h32(tau_vec32) = sqrt(.5)*(randn(Np32,1) + j*randn(Np32,1)); % complex channel path gains % Channel 33 Np33 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec33 = round(L*Ns*rand(Np33,1)); % The location of path delays at the upsampling res. tau_vec33 = tau_vec33 - min(tau_vec33)+1; % Normalized path delays h33 = zeros(L*Ns+1,1); % initialization of channel h33(tau_vec33) = sqrt(.5)*(randn(Np33,1) + j*randn(Np33,1)); % complex channel path gains % Channel 34 Np34 = round(25*rand(1)+1) % Random number of paths for location 11 tau_vec34 = round(L*Ns*rand(Np34,1)); % The location of path delays at the upsampling res. tau_vec34 = tau_vec34 - min(tau_vec34)+1; % Normalized path delays h34 = zeros(L*Ns+1,1); % initialization of channel h34(tau_vec34) = sqrt(.5)*(randn(Np34,1) + j*randn(Np34,1)); % complex channel path gains %-------------- Ni = 100; % Number of iterations for Pe calculations results = zeros(Ni,1); % Initialize correct/incorrect decision vector per each SNR SNRtotals = zeros(ns,1); % Initialize Percentage location error per SNR value for i1 = 1:ns % Loop through all SNR values i1 % Display SNR value currently in play snr_temp = snr(i1); rho = 10^(snr_temp/10); % Convert SNR to dB b = sign(randn(Ni,1)); % Transmitted bit sequence s_F = sqrt(rho)*q; % Reader interrogation signal signal %Inventoried Signals % Location 11 x11 = conv(h11,s_F); % The transmitted signal after passing through the channel Nt = length(x11); % Length of received signal vector n11 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x11 = x11 + n11; % Add AWGN noise onto received signal % Location 12 x12 = conv(h12,s_F); % The transmitted signal after passing through the channel Nt = length(x12); % Length of received signal vector n12 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x12 = x12 + n12; % Add AWGN noise onto received signal % Location 13 x13 = conv(h13,s_F); % The transmitted signal after passing through the channel Nt = length(x13); % Length of received signal vector n13 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x13 = x13 + n13; % Add AWGN noise onto received signal % Location 14 x14 = conv(h14,s_F); % The transmitted signal after passing through the channel Nt = length(x14); % Length of received signal vector n14 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x14 = x14 + n14; % Add AWGN noise onto received signal % Location 21 x21 = conv(h21,s_F); % The transmitted signal after passing through the channel Nt = length(x21); % Length of received signal vector n21 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x21 = x21 + n21; % Add AWGN noise onto received signal % Location 22 x22 = conv(h22,s_F); % The transmitted signal after passing through the channel Nt = length(x22); % Length of received signal vector n22 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x22 = x22 + n22; % Add AWGN noise onto received signal % Location 23 x23 = conv(h23,s_F); % The transmitted signal after passing through the channel Nt = length(x23); % Length of received signal vector n23 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x23 = x23 + n23; % Add AWGN noise onto received signal % Location 24 x24 = conv(h24,s_F); % The transmitted signal after passing through the channel Nt = length(x24); % Length of received signal vector n24 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x24 = x24 + n24; % Add AWGN noise onto received signal % Location 31 x31 = conv(h31,s_F); % The transmitted signal after passing through the channel Nt = length(x31); % Length of received signal vector n31 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x31 = x31 + n31; % Add AWGN noise onto received signal % Location 32 x32 = conv(h32,s_F); % The transmitted signal after passing through the channel Nt = length(x32); % Length of received signal vector n32 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x32 = x32 + n32; % Add AWGN noise onto received signal % Location 33 x33 = conv(h11,s_F); % The transmitted signal after passing through the channel Nt = length(x33); % Length of received signal vector n33 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x33 = x33 + n33; % Add AWGN noise onto received signal % Location 34 x34 = conv(h34,s_F); % The transmitted signal after passing through the channel Nt = length(x34); % Length of received signal vector n34 = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x34 = x34 + n34; % Add AWGN noise onto received signal for i2 = 1:Ni % Repeat transmission simulation Ni times % Randomly select a tag location in the inventoried grid randloc = ceil(12*rand(1)); if randloc == 1 h = h11; end if randloc == 2 h = h12; end if randloc == 3 h = h13; end if randloc == 4 h = h14; end if randloc == 5 h = h21; end if randloc == 6 h = h22; end if randloc == 7 h = h23; end if randloc == 8 h = h24; end if randloc == 9 h = h31; end if randloc == 10 h = h32; end if randloc == 11 h = h33; end if randloc == 12 h = h34; end x_F = conv(h,s_F); % The transmitted signal after passing through the channel Nt = length(x_F); % Length of received signal vector n_F = sqrt(0.5)*(randn(Nt,1) + j*randn(Nt,1)); % Create AWGN noise %x_F = x_F + n_F; % Add AWGN noise onto received signal x_F = x_F/sqrt(sum(x_F.^2)); % normalizes signal values x11 = x11/sqrt(sum(x11.^2)); x12 = x12/sqrt(sum(x12.^2)); x13 = x13/sqrt(sum(x13.^2)); x14 = x14/sqrt(sum(x14.^2)); x21 = x21/sqrt(sum(x21.^2)); x22 = x22/sqrt(sum(x22.^2)); x23 = x23/sqrt(sum(x23.^2)); x24 = x24/sqrt(sum(x24.^2)); x31 = x31/sqrt(sum(x31.^2)); x32 = x32/sqrt(sum(x32.^2)); x33 = x33/sqrt(sum(x33.^2)); x34 = x34/sqrt(sum(x34.^2)); m11 = max(real(conv(x_F,fliplr(x11)))); % calculates max of matched filter m12 = max(real(conv(x_F,fliplr(x12)))); % calculates max of matched filter m13 = max(real(conv(x_F,fliplr(x13)))); % calculates max of matched filter m14 = max(real(conv(x_F,fliplr(x14)))); % calculates max of matched filter m21 = max(real(conv(x_F,fliplr(x21)))); % calculates max of matched filter m22 = max(real(conv(x_F,fliplr(x22)))); % calculates max of matched filter m23 = max(real(conv(x_F,fliplr(x23)))); % calculates max of matched filter m24 = max(real(conv(x_F,fliplr(x24)))); % calculates max of matched filter m31 = max(real(conv(x_F,fliplr(x31)))); % calculates max of matched filter m32 = max(real(conv(x_F,fliplr(x32)))); % calculates max of matched filter m33 = max(real(conv(x_F,fliplr(x33)))); % calculates max of matched filter m34 = max(real(conv(x_F,fliplr(x34)))); % calculates max of matched filter filtervals = [m11,m12,m13,m14,m21,m22,m23,m24,m31,m32,m33,m34]; % Vector of filter outputs [maxvalue,index] = max(filtervals); % Location of max filter output if index == randloc; % Assign 0 if correct location chosen, 1 if not results(i2) = 0; else results(i2) = 1; end end SNRtotals(i1) = (sum(results)/Ni)*100 % Calculate percentage of incorrect location choices per SNR end figure semilogy(snr,SNRtotals) ylabel('Probability of Location Decision Error') xlabel('SNR (dB)') title('RFID Location Decision via Wireless Signatures Error Simulation')
Reply by ●July 24, 20082008-07-24
Perhaps I'm reading your code incorrectly, but based on what I'm seeing, you're not actually implementing spread spectrum of any sort. You seem to simply be taking a random vector (q), and sending it through 12 channels, each with different multipath components. Since each of your signals contains the exact same vector q, you'll simply get a bunch of correlation peaks at each of the multipath delays. Mark
Reply by ●July 24, 20082008-07-24
>Perhaps I'm reading your code incorrectly, but based on what I'm seeing, >you're not actually implementing spread spectrum of any sort. You seemto>simply be taking a random vector (q), and sending it through 12channels,>each with different multipath components. Since each of your signals >contains the exact same vector q, you'll simply get a bunch ofcorrelation>peaks at each of the multipath delays. > >Mark >So are you saying that my problem is the signal (q) rather than the filtering and detection process? That part of the code that creates the signal itself was written by someone else so I just took it for what it was. So would changing the signal solve my problems then or is there more going on? Once again I really do appreciate everyone's help, this is quite a forum I have to say.
Reply by ●July 28, 20082008-07-28
>So are you saying that my problem is the signal (q) rather than the >filtering and detection process? That part of the code that creates the >signal itself was written by someone else so I just took it for what it >was. So would changing the signal solve my problems then or is theremore>going on? Once again I really do appreciate everyone's help, this isquite>a forum I have to say.Sorry it took so long to reply, but I've been vacationing. In answer to your question, sort of. What your "signal" seems to be is nothing more than a random vector of bits (BPSK), sent through 12 different channels, each with unique multipath content (gain and delay, as well as number of paths). What you can do with this signal set is determine the multipath channel impulse response for each case, which will show up as correlation peaks for each component (with corresponding amplitude, phase and delay). If we assume that your "signal" code is correct for what you're really trying to do, then your detection problem is really a channel identification problem. What you would need to do here would be to take your convolution output, pick the largest peaks, then compare them to the known multipath channels to see which it is. Picking the largest peaks, btw, is really an arbitrary exercise. One method would be to find the mean value of the convolution output (squared, for a sample variance comparison), then choose any peaks that have a greater squared value than some arbitrary level above the mean, say "choose all peaks that are 6 dB above the mean." In radar, this would be a noise-riding threshold detector. The list of peaks (amplitude, phase and delay) would then serve as your impulse response vector, and you would compare to each of the templates (known channels) to determine which is closest to the signal received. Spread spectrum, however, is a bit more than just this. The simplest to implement is direct sequence spread spectrum (DSSS), and is what most PCS (cellular) phones in the US use (CDMA, actually). DSSS functions by taking your random vector (called a spreading code), q, and using that to "spread" the data into a wider bandwidth than the data by itself. This is done simply by XORing each bit of data (assume BPSK for simplicity) with q, and transmitting at length(q) times the bit-rate. In your example, the transmitted bandwidth would then be spread by a factor of 128. A string of 5 data bits such as [1, -1, -1, 1, 1] would be transmitted as [q, -q, -q, q, q] where the latter is 5*length(q) bits. Switching back to your desired application, what you might be trying to do is to identify which of N (I suppose N = 12 in this example) spreading sequences is being used in a received signal. What you would then do would be generate N different spreading codes q_i, i = 1, ..., N. Normally these codes are taken from a set of known orthogonal (uncorrelated) sequences, such as Gold or Kasami sequences since this minimizes cross-correlations between separate users, thus making detection more robust. You might perform the convolution in the same manner, trying each of your N codes (q_i), while looking for the code that provides correlation peaks every 128 samples. Multipath will blur your outputs, btw, particularly if you aren't using orthogonal codes. Mark
