Tim Wescott wrote:
> On Sat, 23 Feb 2008 04:06:57 -0800, bhargavajs.07 wrote:
>
>> Hi all
>> I'm working on GMSK modem simulation in MATLAB. I'm almost done with
>>
>> the modulation part albeit i'm not sure if the model I used is correct
>>
>> adding to my trouble is the ambiguity with the Power Spectral Density
>>
>> of the GMSK modulated signal (also the MSK modulated signal). Here I'm
>>
>> writing the algorithm I followed to do the simulation:
>>
>> bitrate=4000bps BT=0.3
>>
>> 1.calculated the impulse response of the GaussianLPF with the formula
>>
>> g(t)= convolution(gaussian pulse , rectangular pulse of width T) as
>>
>> specified in the GSM 5.04 version 8.1.2 Release 1999. In fact I directly
>> used the impulse response specified in the
>>
>> following link
>>>> http://www.emc.york.ac.uk/reports/linkpcp/appD.pdf
>> 2. then I sampled the incoming binary data stream exactly in the
>>
>> middle of the bit period, then over sampled the stream with an over
>>
>> sampling rate of 41
>> (here in a bit period first 20 samples will be 0 the 21st sample will
>>
>> be +/-1 depending on the data thenext 20 will be 0 again)
>>
>> 3. I convolved these samples with the Impulse response I calculated
>>
>> earlier.
>>
>> 4.then integrated the resulting signal to get 'y'
>>
>> 5.then calculated the I and Q channel streams as I=cos(y),Q=sin(y)
>>
>> then modulated it with a I/Q modulator with Fc(carrier frequency=5000)
>>
>> then I got the modulated signal as shown in the link below for the
>>
>> data:[1 1 1 1 1 0 0 0 0 0]
>>>> http://aycu19.webshots.com/image/45298/2003374167469833630_rs.jpg
>> Does a GMSK modulated signal look like that??? It shows constant
>> amplitude and frequency deviations the Ichannel vs Q channel (baseband)
>> plot is as shown in the following
>>
>> link
>>>> http://aycu15.webshots.com/image/45574/2003383907486945009_rs.jpg
>> It is a circle with radius 1 as is expected of a GMSK signal
>>
>>
>> I calculated its PSD of the above signal using FFT and scaled its
>>
>> frequency axis to 0 to Fs/2 (here Fs=16400Hz) It is in the zoomed in
>> figure(for viewability) shown in following link
>> http://aycu20.webshots.com/image/46459/2003353100129112203_rs.jpg
>>
>> My doubt here is:
>> In many text books and technical papers the PSD of GMSK signal shows
>>
>> a peak at the carrier frequency but in my signal's PSD shows a dip
>>
>> there I'm not getting where I went wrong.
>>
>> PLEASE GUIDE ME WHERE I WENT WRONG AND PLEASE TELL ME THE CORRECT WAY
>>
>> OUT OF THIS PROBLEM
>
> Your plots are entirely as expected. GMSK, like MSK, is a form of
> frequency shifting: the MSK part specifies that the phase will be
> continuous and that the frequency shift will be exactly 1/2 the bit rate;
> the 'G' means that the frequency shift will be filtered with a Gaussian
> profile before modulating the carrier. So sending a bunch of ones then a
> bunch of zeros will give you one tone, then another.
>
> Your PSD calculation is a consequence of your bit pattern: in using a
> burst of ones followed by a burst of zeros you have failed to meet the
> requirements for calculating a PSD, which is to use a random bit pattern
> as John mentioned. Using a random bit pattern, you will find, requires
> that you use a long interval to make sure that there is little
> correlation between the bits. Even so, a PSD that's calculated
> experimentally will be rougher than one that's calculated from first
> principals (and no, I can't tell you how to do this from first
> principals: I haven't taken the time to do it myself, and I haven't found
> a paper that discusses it).
>
I just received the following email, I thought it deserved a wider audience:
Hello Tim,
My name is Mick Owens. I stumbled upon your response to the poor
soul that is trying to recreate the GMSK figures from that paper from
the University of Hull. I too have gone down that path and after
several missteps, I finally got it.
Then I took it one step further and plotted the PSD of the
modulated GMSK signal. It came out text-book perfect, so I thought I'd
share it with you. The attached Matlab file is fairly well commented, so
you shouldn't have much trouble with it. The file will plot each of the
University of Hull figures, although they are all commented out in the
current version to optimize the performance of the spectral analysis.
I could not determine the email address of the person you were
responding to, so I would appreciate if you could send him the attached
file.
Thanks, Mick
% 3/19/2008
% M. Owens
%
% This routine implements the GMSK encoding algorithm described in the
University of Hull
% paper called "Appendix D - Digital Modulation and GMSK"
%
% Create the same repeating string of ones and zeros that is defined in
the Hull paper
%data_bits = [1,1,-1,1,1,-1,-1,1,-1,1,-1,-1,1,1,-1,1,1,-1,-1,1,-1,1,-1,-1];
%
% SPECTRAL ANALYSIS ONLY: Create a very long random string for spectral
analysis
size_of_data = 1000;
data_bits = zeros(1,size_of_data);
for n = 1 : size_of_data
if rand < 0.5
data_bits(n) = 1;
else
data_bits(n) = -1;
end
end
% END SPECTRAL ANALYSIS ONLY
%
% Define the modulation parameters
BT = 0.5; % This Bandwidth-Time product affects the Inter-Symbol
Interference (ISI)
R = 265; % Bit Rate, bps
fc = 5000; % Carrier frequency, Hz
time_steps = 250; % This sets the sampling rate of the modulated signal
%
% Pre-calculate several variables to save time in later iterations:
T = 1/R; % Bit time, seconds
B = BT/T; % 3 dB Gaussian filter bandwidth
C1 = 2*pi*BT/(T*sqrt(log(2)));
C2 = 1/(2*T);
C3 = T/2;
%
% Define the time vector for evaluating the Q-function. This performs the
% truncation described in the Hull paper.
t = -T : T/time_steps : T; % Note the use of "negative time"
%
% Break the Q-function into two parts, plus and minus
tau_minus = C1*(t - C3); % Avoid the use of a for loop here using t
vector. Saves time in Matlab.
Q_minus = 0.5*erfc(tau_minus/sqrt(2)); % Equivalent to Q-function
defined in Hull paper
tau_plus = C1*(t + C3);
Q_plus = 0.5*erfc(tau_plus/sqrt(2));
%
% Develop a normalization factor to keep the peak amplitude at a
constant '1'
tau_nml_minus = C1*(0 - C3);
Q_nml_minus = 0.5*erfc(tau_nml_minus/sqrt(2));
tau_nml_plus = C1*(0 + C3);
Q_nml_plus = 0.5*erfc(tau_nml_plus/sqrt(2));
Q_nml = (1/(2*T))*(Q_nml_minus - Q_nml_plus);
%
% Compute the Gaussian response to a single positive impulse. This
creates the
% template that will be used to string together a series of responses later.
g = (1/(2*T))*(Q_minus - Q_plus)/Q_nml;
%
% The following plot recreates the filter response shown in the Hull paper
% plot(g)
%
% Create the b(t) array to hold the series of Gaussian responses
b_t = zeros(1, (length(data_bits)+1)*time_steps+1);
%
% Add the individual impulse response templates to b(t) at the time
locations and polarities
% defined by the data_bits
for n = 1 : length(data_bits)
temp_array = zeros(1, (n-1)*time_steps); % This adds leading zeros
before the template
temp_array2 = [temp_array data_bits(n)*g]; % Put a pos or neg
template into the array
temp_array3 = zeros(1, length(b_t) - length(temp_array2)); % Create
trailing zeros
temp_array4 = [temp_array2 temp_array3];
%
% Recreate the individual shaped pulses shown in the Hull paper:
% hold on
% plot(temp_array4)
%
b_t = b_t + temp_array4;
end
%
% The following plot recreates the function b(t) shown in the Hull paper
% plot(b_t)
%
% Integrate b(t) to create c(t) as described in Hull
for n = 2 : length(b_t)
b_t(n) = b_t(n - 1) + b_t(n); % This is a very simple square-rule
integration
end
c_t = b_t; % Transfer the data to keep the same variable name as Hull
%
% Recreate the function c(t) shown in the Hull paper
% plot(c_t)
%
% The response of the filter to a single '1' in the data stream should
create a phase change
% of pi/2. But the I and Q functions use cos and sin to calculate their
arguments in radians,
% so a scale factor is necessary. To determine the scale factor, sf, use
a single impulse
% response and take the last value in c_t.
temp_t = g; % use a single impulse response
% Integrate the single bit response
for n = 2 : length(temp_t)
temp_t(n) = temp_t(n - 1) + temp_t(n); % Square-rule integration
end
% Generate the scale factor, sf, using the last element in temp_t
sf = pi/(2*temp_t(length(temp_t)));
I = cos(sf.*c_t);
Q = sin(sf.*c_t);
%
% The following plots recreate the figures of I & Q shown in the Hull paper
% plot(I)
% plot(Q)
%
I = I * pi/2; % Convert the I & Q signals into radians
Q = Q * pi/2;
%
% Create the time vector, t, used to create the modulated signal
cycle_time = 1/fc; % Length in seconds
cycles = length(data_bits) * fc/R; % # bits x cycles/bit = cycles
end_time = cycles * cycle_time;
samples_sec = 1/(end_time/length(c_t));
t = 0 : end_time/length(c_t) : end_time - end_time/length(c_t);
%
% Create the modulated signal using the following trig identity:
% cos(x+ph) = cos(x)cos(ph) - sin(x)sin(ph)
%
% where: x = 2*pi*fc*t
% cos(ph) = phase contribution of I (i.e., cos(c_t))
% sin(ph) = phase contribution of Q (i.e., sin(c_t))
% sin(x+ph) = transmitted signal, including I & Q phase modulation
%
mod_signal = cos(2*pi*fc*t) .* I - sin(2*pi*fc*t) .* Q;
% Note that the '.*' operator performs an element-by-element
multiplication in Matlab
% plot(mod_signal);
%
% For spectral analysis, use Matlab's built-in FFT function. Start by
% applying a Hamming window to the signal samples in 'mod_signal'.
size_of_mod_sig = length(mod_signal);
% Apply the Hamming window to 'mod_signal'. This smooths the amplitude
modulation that would
% otherwise affect the FFT results if the sample window was simply
started and stopped.
for n = 1 : 1 : size_of_mod_sig
mod_signal(n) = mod_signal(n) * (0.53836 - (0.46164 *
cos((2*pi*n)/(size_of_mod_sig - 1))));
end
%
% Now take the FFT
FFT_points = 2^nextpow2(length(mod_signal));
Y = abs(fft(mod_signal,FFT_points)).^2; % abs gives magnitude, .^2
converts voltage to power: V2/R
Y = 10*log10(Y/.001) - 100; % Convert to dBm and normalize the Y-scale
to 0 dBm
%
% Prep for graphing:
freq_per_bin = samples_sec/FFT_points; % frequency per FFT bin
bins = FFT_points/2; % number of bins before the foldover at midpoint
(Nyquist freq)
hf = freq_per_bin * bins; % highest frequency on the x-axis
%
% Just interested in the GMSK spectrum around the center frequency, so zoom
% in on that data:
center_freq_bin = fc/freq_per_bin;
graph_width = 5 * R / freq_per_bin;
high_bin = round(center_freq_bin + graph_width / 2);
low_bin = round(center_freq_bin - graph_width / 2);
graph_bins = high_bin - low_bin;
f = freq_per_bin*(low_bin:high_bin);
plot(f,Y(low_bin:high_bin))
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density (dB)');
--
Tim Wescott
Wescott Design Services
http://www.wescottdesign.com
Do you need to implement control loops in software?
"Applied Control Theory for Embedded Systems" gives you just what it says.
See details at http://www.wescottdesign.com/actfes/actfes.html