FSK Demod using Quadrature Receiver with matched filters
in [DSP]
|
Prev: whitening filter
Next: GMSK Simulation
From: GSB on 5 Nov 2005 18:55 I created a Matlab M-file for a FSK demodulator using a Quadrature Receiver with 4 matching filters. The code is found below. I need someone to verify that this code works or to point out any problems. Also, if there is anything I can do to make it stronger and better, then I would like to hear about it. My current problem is that awgn(in_signal, 18, 'measured') of 18 is the lowest SNR I can use before I get bit error rates. Is there a way to increase this codes ability to remove noise while decoding the signal? I have tried lowpass and bandpass filters. But this only produces more bit error rates. Also, how can I create an automatic gain control (AGC)? % Keep function result = test % test4_quadrature_receiver() do_loop = 1; iter = 1; while do_loop == 1 [result1 result2 result3] = test1(); result3 = 1; % Do not test this if result1 ~= 20 || result2 ~= 20 || result3 ~=1 disp('ERROR Found!!!'); do_loop = 0; else disp(sprintf('PASSED -- Iteration #%d', iter)); end iter = iter + 1; % do_loop = 0; % iterate only 1 time end % Keep function [out1 out2 out3] = test1 clear all; set_rand_seed(); % Create Signal symbols1 = [1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1]; pattern1 = [2 2 2 1 2 1 1 2 2 1]; num_patterns1 = 20; for i = 1:num_patterns1 symbols1 = [symbols1, pattern1]; end len_symbols1 = length(symbols1); symbols2 = [1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1]; pattern2 = [1 2 3 4 2 3 1 4 3 4 2 1 2 1 4 3]; num_patterns2 = 20; for i = 1:num_patterns2 symbols2 = [symbols2, pattern2]; end pattern3 = [1 2 3 4 4 3 2 1 1 2 3 4 4 3 2 1]; symbols2 = [symbols2 pattern3 [1 2 3 4 1 1 2 2 1 3 4 4 2 1 2 2 1 3 3 4 2 1 2 3 4 1 1 2 2 1 3 4 4 2 1 2 2 1 3 3 4 2 1 2 3 4 1 1 2 2 1 3 4 4 2 1 2 2 1 3 3 4 2]]; min_num_samples = 18; % should be a little lower than samples_per_bit Fc = 31500; Fs = 4*Fc; bps = floor(get_rand_num(7000, 1000)); frequencies = [Fc - 800, Fc + 800, Fc - 2400, Fc + 2400]; samples_per_bit = 26; % ceil(Fs / 4800) do_loop = 0; while do_loop == 1 samples_per_bit = floor(3*min(frequencies) / bps); t = ceil(min_num_samples * 1.10); % Needs to be 10% larger since real minimum might be smaller if samples_per_bit < 26 disp(sprintf('ERROR: samples_per_bit is %d but must be >= %d', samples_per_bit, t)); disp(sprintf(' : Either increase Fc and/or decrease bps. Must be (3*Fc) / 13 >= bps.')); else do_loop = 0; end end % Symbols to signal num_samples1 = 50; t1 = symbols_to_signal(symbols1, frequencies, samples_per_bit, Fs); t2 = symbols_to_signal(symbols2, frequencies, samples_per_bit, Fs); % min_num_samples = min(samples_per_bit, num_samples1); symbols1 = []; symbols2 = []; in_signal = [t1 t2]; % Noise in_signal = awgn(in_signal, 18, 'measured'); % Filter init_test5_quadrature_receiver init_quadrature_receiver_per_bit(); init_remove_small_transitions_in_symbols(); init_symbols_to_bits(); len = length(in_signal); num_bit_per_iter = floor(get_rand_num(20, 2000)); disp(sprintf('bps: %d, samples_per_bit: %d, min_num_samples: %d, num_bit_per_iter: %d', bps, samples_per_bit, min_num_samples, num_bit_per_iter)); out_symbols = []; for i = 1:num_bit_per_iter+1:len - num_bit_per_iter t = test5_quadrature_receiver(in_signal(i:i+num_bit_per_iter), frequencies, min_num_samples, Fs); out_symbols = [out_symbols t]; end if (i + num_bit_per_iter) < len t = test5_quadrature_receiver(in_signal(i+num_bit_per_iter+1:len), frequencies, min_num_samples, Fs); out_symbols = [out_symbols t]; end % BERT - How many patterns found? out1 = total_num_pattern(out_symbols, pattern1); out2 = total_num_pattern(out_symbols, pattern2); out3 = total_num_pattern(out_symbols, pattern3); disp(sprintf('Found %d/%d, %d/%d, %d/%d', out1, num_patterns1, out2, num_patterns2, out3, 1)); % Keep function out_bits = init_test5_quadrature_receiver() global g1_in_signals g1_decoded_symbols g1_clean_symbols g1_out_bits g1_in_signals = []; g1_decoded_symbols = []; g1_clean_symbols = []; g1_out_bits = []; % Keep function out_bits = test5_quadrature_receiver(in_signal, frequencies, min_samples_per_bit, Fs) global g1_in_signals g1_decoded_symbols g1_clean_symbols g1_out_bits min_num_data_to_process = 10; out_bits = []; % Filters % Demodulate g1_in_signals = [g1_in_signals in_signal]; if length(g1_in_signals) < min_num_data_to_process % disp('d1'); return; end decoded_symbols = signals_to_symbols_via_quadrature_receiver(frequencies, g1_in_signals, min_samples_per_bit, Fs); g1_in_signals = []; % clean signals g1_decoded_symbols = [g1_decoded_symbols decoded_symbols]; if length(g1_decoded_symbols) < min_num_data_to_process % disp('d2'); return; end clean_symbols = remove_small_transitions_in_symbols(g1_decoded_symbols, 5); g1_decoded_symbols = []; % symbols to bits g1_clean_symbols = [g1_clean_symbols clean_symbols]; if length(g1_clean_symbols) < min_num_data_to_process % disp('d3'); return; end num_alt_bits = 10; out_bits = symbols_to_bits(g1_clean_symbols, num_alt_bits); g1_clean_symbols = []; % Keep function result = plot_spectrum(signal, Fs) result = []; plot(Fs*(0:1/(length(signal)-1):1), abs(fft(signal))); % Keep function result = set_rand_seed() rand('seed', sum(100*clock())); result = 0; % Keep function result = get_rand_num(min, max) result = ((max - min) * rand(1)) + min; %Keep function result = init_quadrature_receiver_per_bit() global g_quad_last_bits; g_quad_last_bits = []; % Keep % DESCRIPTION: This function uses Quadrature Receiver techniques to % determine which of the reference frequencies the signal matches closest. % The entire signal is sent which contains multiple symbols. It handles % any number of reference frequencies. % FUNCTION: quadrature_receiver(referencies, signals, Fs) % where: referencies - All reference frequencies, e.g. [f1 f2 f3 f4]. % signals - Signal data, e.g. sin wave % Fs - Number of samples per second. % SYNTAX: quadrature_receiver([100 200 300 400], sin([0:0.02:2*pi]), 5000) % result = One of the indexes of frequencies corrsponding to closest % frequency the signal matches. % function result = signals_to_symbols_via_quadrature_receiver(frequencies, signals, min_samples_per_bit, Fs) global g_quad_last_bits; min_samples_per_bit = min_samples_per_bit; % 10 result = []; signals = [g_quad_last_bits signals]; len_signals = length(signals); for i = 1:1:(len_signals - min_samples_per_bit + 1) a = quadrature_receiver_per_bit(frequencies, signals(i:i+min_samples_per_bit-1), Fs); result = [result, a]; end g_quad_last_bits = signals(len_signals - min_samples_per_bit + 2:len_signals); % Keep % DESCRIPTION: This function uses Quadrature Receiver techniques to % determine which of the reference frequencies the signal matches closest. % The entire signal is compared to all of the reference frequencies. It % handles any number of reference frequencies. % FUNCTION: quadrature_receiver_per_bit(referencies, signal, Fs) % where: referencies - All reference frequencies, e.g. [f1 f2 f3 f4]. % signal - Signal data, e.g. sin wave % Fs - Number of samples per second. % SYNTAX: quadrature_receiver_per_bit([100 200 300 400], sin([0:0.02:2*pi]), 5000) % result = One of the indexes of frequencies corrsponding to closest % frequency the signal matches. % function result = quadrature_receiver_per_bit(frequencies, signal, Fs) z_values = []; len = length(signal); for i = 1:length(frequencies) ref_cos(i,:) = get_samples_from_freq('cos', frequencies(i), len, Fs); ref_sin(i,:) = get_samples_from_freq('sin', frequencies(i), len, Fs); z1(i) = sum(sqrt(2/len) * (ref_cos(i, 1:len)) .* signal(1:len)); z2(i) = sum(sqrt(2/len) * (ref_sin(i, 1:len)) .* signal(1:len)); z_values = [z_values, z1(i)^2 + z2(i)^2]; end result = find(z_values == max(z_values)); result = result(1); % disp(sprintf('min_len1: %d, min_len2: %d, z1: %d, z2: %d, z3: %d, z4: %d, result: %d', min_len1, min_len2, z1, z2, z3, z4, result)); % Keep function result = total_num_pattern(symbols, pattern) len_pattern = length(pattern); result = 0; for i = 1:length(symbols) - len_pattern + 1 if symbols(i:i+len_pattern-1) == pattern result = result + 1; i = i + len_pattern - 1; end end %Keep function result = init_remove_small_transitions_in_symbols() global num_no_transitions last_signal num_no_transitions = 0; last_signal = []; %Keep function out_signal = remove_small_transitions_in_symbols(in_signal, min_transition_bit) global num_no_transitions last_signal; out_signal = []; in_signal = [last_signal in_signal]; for i = 2:length(in_signal) % if no transition, skip if in_signal(i) == in_signal(i-1) num_no_transitions = num_no_transitions + 1; continue; end % if time is long enough, then get transition if num_no_transitions >= min_transition_bit len = length(out_signal); out_signal(len+1:len+num_no_transitions+1) = in_signal(i - 1); num_no_transitions = 0; end end last_signal = in_signal(length(in_signal)); %Keep function result = get_samples_from_freq(method, freq, num_samples, Fs) p = freq * ((0.2 / pi) * (100/Fs)); angle = [p*ones(1,num_samples)]; phase = cumsum(angle) - angle(1); if strcmp(method, 'sin') == 1 result = sin(phase); elseif strcmp(method, 'cos') == 1 result = cos(phase); else disp('Error: method var is undefined in get_1hz_from_freq(), exiting'); exit; end % Keep function result = symbols_to_signal(symbols, frequencies, samples_per_symbol, Fs) % Usage: a = symbols_to_signal([1 2], [10 20], 36000, 150000); % plot(a); % Notes: % - Symbols is a 1-n array with values 1..n for i = 1:length(symbols) p = frequencies(symbols(i))* ((0.2 / pi) * (100/Fs)); t = [p*ones(1,samples_per_symbol)]; if i == 1 angle = t; else angle = [angle, t]; end end phase = cumsum(angle) - angle(1); result = sin(phase); % Keep function init_symbols_to_bits() global last_bit cur_bit cur_num_alt_bits bit_time first_bit_time sum_bit_time gindex cur_bit_time last_bit_time hold_out_bits is_initialized; is_initialized = 1; bit_time = 999999; first_bit_time = 123456; sum_bit_time = 0; gindex = 0; cur_bit_time = 1; cur_num_alt_bits = 0; cur_bit = 0; hold_out_bits = []; % Keep function out_bits = symbols_to_bits(in_bits, num_alt_bits) global last_bit cur_bit cur_num_alt_bits bit_time first_bit_time sum_bit_time gindex cur_bit_time last_bit_time hold_out_bits is_initialized; % Test if init_symbols_to_bits() called. Does not work, but fix it. Test % if bit_time was defined and set to 999999. if is_initialized ~= 1 disp('Need to call init_symbols_to_bits(), exiting'); return; end out_bits = []; for i = 1:length(in_bits) gindex = gindex + 1; last_bit = cur_bit; cur_bit = in_bits(i); % Sample Bit t = (gindex - cur_bit_time) / round(bit_time / 2); if (floor(t) == t) && (floor(t/2) ~= (t/2)) % disp(sprintf('bit sample: %d', cur_bit)); out_bits = [out_bits, cur_bit]; end % No transitions if cur_bit == last_bit continue; end % Do not let gindex become too large if cur_bit_time > 10000 cur_bit_time = cur_bit_time - 8000; gindex = gindex - 8000; end % Update bit times last_bit_time = cur_bit_time; cur_bit_time = gindex; % If transition is same bit time. t1 = abs(cur_bit_time - last_bit_time - first_bit_time); t2 = first_bit_time * 0.20; if t1 <= t2 || (first_bit_time == 123456) % First bit in alternating bits if cur_num_alt_bits == 0 first_bit_time = cur_bit_time - last_bit_time; end sum_bit_time = sum_bit_time + (cur_bit_time - last_bit_time); cur_num_alt_bits = cur_num_alt_bits + 1; % disp(sprintf('t1: %d, t2: %d, cur_bit_time: %d, last_bit_time: %d, sum_bit_time: %d, first_bit_time: %d, cur_num_alt_bits: %d', t1, t2, cur_bit_time, last_bit_time, sum_bit_time, first_bit_time, cur_num_alt_bits)); else cur_num_alt_bits = 0; sum_bit_time = 0; first_bit_time = 123456; % disp(sprintf('t1: %d, t2: %d, cur_bit_time: %d, last_bit_time: %d, sum_bit_time: %d, first_bit_time: %d, cur_num_alt_bits: %d', t1, t2, cur_bit_time, last_bit_time, sum_bit_time, first_bit_time, cur_num_alt_bits)); end % if got 10 bit transitions if cur_num_alt_bits >= num_alt_bits % disp('************Found 10 alternating bits'); bit_time = round(sum_bit_time / cur_num_alt_bits); cur_num_alt_bits = 0; first_bit_time = 123456; sum_bit_time = 0; end end % disp(sprintf('gindex: %d', gindex)); % disp(sprintf('bit_time: %d', bit_time));
From: Tim Wescott on 5 Nov 2005 19:38 GSB wrote: > I created a Matlab M-file for a FSK demodulator using a Quadrature > Receiver with 4 matching filters. The code is found below. I need > someone to verify that this code works or to point out any problems. > Also, if there is anything I can do to make it stronger and better, > then I would like to hear about it. > > My current problem is that awgn(in_signal, 18, 'measured') of 18 is > the lowest SNR I can use before I get bit error rates. Is there a way > to increase this codes ability to remove noise while decoding the > signal? I have tried lowpass and bandpass filters. But this only > produces more bit error rates. > > > Also, how can I create an automatic gain control (AGC)? > <lots - o - code snipped> That's a lot of code for someone to wade through for free -- I suspect it'd take a good half hour to two hours to grok what you have there, and I couldn't begin to say how long to find any problems. If no one takes this on I suggest that you reduce the code down to a block diagram or an algorithmic description, and try those out on us. In the process you may see your problem. In general you can increase a detector's performance by using matched filters, assuming that you have truly orthogonal signals. If the FSK is phase continuous and the phase changes are well behaved enough you may be able to treat it as a trellis-coded modulation scheme. If your SNR is really that bad when you start seeing problems then I would suspect you have a plain old bug, however. Have you generated an eye diagram to see what the result looks like? You should have nice wide open eyes with noiseless data; if you don't then you have bugs or intersymbol interference. Jaggy or jumpy lines through the eyes would indicate logic bugs; smooth ones could either be bugs or intersymbol interference. You create an automatic gain control by measuring the magnitude of the signal at some point in the process and applying a varying gain (in the algorithmic case you just set a multiplier). If you're doing this as a batch process you just need to normalize your data set to your desired power level before demodulation. If you are working with a real receiver that has a hardware gain change then you need to measure some representative signal in your processing chain and adjust the gain to keep your desired signal constant. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com
From: Bevan Weiss on 5 Nov 2005 22:41 GSB wrote: > I created a Matlab M-file for a FSK demodulator using a Quadrature > Receiver with 4 matching filters. The code is found below. I need > someone to verify that this code works or to point out any problems. > Also, if there is anything I can do to make it stronger and better, > then I would like to hear about it. > > My current problem is that awgn(in_signal, 18, 'measured') of 18 is > the lowest SNR I can use before I get bit error rates. Is there a way > to increase this codes ability to remove noise while decoding the > signal? I have tried lowpass and bandpass filters. But this only > produces more bit error rates. > > > Also, how can I create an automatic gain control (AGC)? I'm definitely going to agree with the others regarding code size... try and keep it to just the relevant points. Also, I'm sure you can condense up your [ 1 2 1 2 1 2...] array declaration heaps by using repeated subarrays. Creating a linearly increasing/decreasing array is also easy using colon notation. ie [1:100] creates 100 values from 1 to 100. You should perhaps use the communications toolbox for these kind of simulations, much easier if you then want help with them. Not quite so much code to trawl through. If you don't have access to it, then blocking it up would probably be a good idea all the same. You should get bit error rates regardless of the signal to noise ratio. Just with a lower snr you will get more errors. So you'd need to provide some quantitative value of errors you're experiencing at 18dB SNR, might be as expected. As has already been stated, you'll want to use matched filters. Though if you're doing non-coherent demodulation then it's a different situation..
From: GSB on 7 Nov 2005 17:18 I tried this stuff using Simulink, but was told to use M-file for more control and accuracy. I had thought asking for help would be a fair trade since I gave the FSK code out for free in exchange for some advice. I suspect there are a few people copying the code for their own use. Essentially, the code works until I add some gaussian noise (AWGN) above 18 dB or so. The code uses a I/Q Quadrature Receiver with Matched Filters based on Sylar's "Digital Communications" book. The stuff works with less noise and at slower speeds. But when I add some noise and use faster baud speeds, I get some errors. For now, I just want to verify I implemented the Matched Filters correctly for non-coherent FM signals. What these matched filters do is sample say 50 bits, at position i, of data and take the highest magnitude and output a symbol which denotes the frequency being detected. Next, I repeat the same process but start the 50 samples at i + 1 position. Some books sample at every 50 bit increments if the samples per bit is 50 bits per symbol. However, my code samples 50 or so bits in increments of 1. Am I doing things correctly here? I am just trying to get a better BER while increasing noise and baud speed. As far as AGC (gain control), how do you get more gain in DSP? I am thinking of measuring the amplitude of each sample/bit or hertz unit of signal and multiplying this portion of the signal by 1/amplitude in order to get the amplitude back to +1 and -1 levels. Thanks,
From: Bevan Weiss on 7 Nov 2005 17:31
GSB wrote: > I tried this stuff using Simulink, but was told to use M-file for more > control and accuracy. I had thought asking for help would be a fair > trade since I gave the FSK code out for free in exchange for some > advice. I suspect there are a few people copying the code for their > own use. Essentially, the code works until I add some gaussian noise > (AWGN) above 18 dB or so. The code uses a I/Q Quadrature Receiver with > Matched Filters based on Sylar's "Digital Communications" book. The > stuff works with less noise and at slower speeds. But when I add some > noise and use faster baud speeds, I get some errors. For now, I just > want to verify I implemented the Matched Filters correctly for > non-coherent FM signals. What these matched filters do is sample say > 50 bits, at position i, of data and take the highest magnitude and > output a symbol which denotes the frequency being detected. Next, I > repeat the same process but start the 50 samples at i + 1 position. > Some books sample at every 50 bit increments if the samples per bit is > 50 bits per symbol. However, my code samples 50 or so bits in > increments of 1. Am I doing things correctly here? I am just trying > to get a better BER while increasing noise and baud speed. > > As far as AGC (gain control), how do you get more gain in DSP? I am > thinking of measuring the amplitude of each sample/bit or hertz unit of > signal and multiplying this portion of the signal by 1/amplitude in > order to get the amplitude back to +1 and -1 levels. > > Thanks, Your matched filter should be an integrate and dump rather than the maximum value (if you're dealing with binary pulse data signals). You may also be running into some issues with threshold effect with the FSK signal if you're not extracting all signal power efficiently (such as not using an integrate and dump matched filter). This occurs when the SNR becomes low, and hence the receiver can no longer synchronise to the frequency changes. Most people on this forum are capable of producing their own FSK receivers, infact many have probably got some in commercial products at this time. So the issue is mostly where the problem lies in yours and exposing just that portion to us so that we don't have to trawl through the rest of it. |