Some fitting problems are intrinsically non-linear and the techniques we have to use to for these are quite different from the linear and polynomial fitting we looked at in the previous exercise. In particular it is important to be able to fit a broad range of non-linear problems. Today we will focus on fitting Gaussian pro- files, which are frequently seen in experimental data in many different contexts, however, the techniques we will investigate can be applied equally well to other nonlinear fitting problems. The formula for a Gaussian profile is
where A is the amplitude of the Gaussian, b is the offset of the profile from x = 0 and c is related to the width of the Gaussian profile.
For today’s exercise we will take some data collected with the University of Tasmania 26m radio telescope. In the KYA320 lectures over the next few weeks we will be studying antennas. One of the properties of real antennas is that they have directional sensitivity – their transmission (or reception) is distributed over a solid angle in a manner which depends on the specific antenna design. This is known as the beam pattern. For a parabolic reflecting antenna (such as the radio telescope) to first order the beam pattern is a Gaussian profile with an angular full-width at half-maximum (FWHM) of
where λ is the wavelength at which the antenna is being operated and D is the antenna diameter.
We will use the collected data to measure the FWHM of the 26m radio telescope and compare our result with the prediction given by equation above. NOTE: the FWHM is related to the Gaussian parameter c, but it isn’t they are not equal. What is the relationship between c and the FWHM? (hint: solve
1. Initial examination of the data
1. On the KYA320 MyLO site you will find two MATLAB files exercise5a.mat and exercise5b.mat, download these files, load them into MATLAB (using load exercise5a.mat). Each of these files contains a scan with the 26m telescope across a “point source”, known as 30 Doradus (a region of ionised gas in the Large Magellanic Cloud). When the telescope is scanned across a point source, the change in signal level observed reflects the angular changes in the telescope beam pattern. The first (data(:,1)) and second (data(:,2)) columns of the data matrix hold the power from one of the radio telescope’s receivers.
2. Have a look at the variables which have been defined and the values that they hold, see if you can work out what they are. For exercise5a.mat the data variable holds the 127050 samples for each of 2 channels. The data was sampled at a frequency of approximately 847 Hz (the exact value is stored in freq), and so 127050 samples corresponds to a time interval of approximately 2.5 minutes.
3. Plot the data for each of the two channels on the same figure. The two sudden changes in the recorded level at the start of the scan are where the receiver “noise diode” (a calibration tone) has been switched on and the addition of 1 dB of attenuation respectively. These measurements are required to properly calibrate the data. Can you see the “source” in these scans (it is quite hard in channel 1)? Do you think that the two channels recorded data at similar frequencies – if not, why not?
4. The raw data from the telescope is sampled rapidly so that any short bursts of radio frequency interference (RFI) can be more easily removed from the data. Do you see any signs for RFI in the data we are analysing and do the aberrant samples have any distinguishing characteristics?
5. The “RFI” in these data will not have any serious impact on our analysis (it is low amplitude and short duration), so we don’t need to worry about excluding it from our analysis. However, it is probably desirable to exclude the noise diode and 1 dB step changes in level from our analysis. One relatively easy way to do this is through the use of either logical arrays or index arrays. A logical array contains an element corresponding to each data point, if the logical array element is 0 then that data point is excluded, if it is 1 then it is included. For example to create a logical array to exclude the noise diode section of the data we can do something like
i n cl u d e 1 = logical ( ones (length ( data ) , 1 ) ) ;
i n cl u d e 1 ( 2 1 7 0 0: 2 5 0 0 0 ) = 0 ;
p lot ( data ( i n cl u d e 1 , 1 ) )
An index array is a array containing a list of the elements from a larger array which you want to include – it can either exclude elements you don’t want, or repeat them with other valid values. To create an index array which performs the same function as the logical array above we want to create an array which contains all the numbers 1 to 127500 once, excluding those between 21700 and 25000. Write some MATLAB code to do this and compare the results with the logical array method.
6. For the purposes of analysis we will sometimes (as in this case) want to perform some sort of averaging on the data to reduce the noise level in an individual data point and to reduce the total number of data points passed to the fitting and other routines. There are many different ways of smoothing the data, one is to simply average up some number of consecutive samples and write out a new data array of the averaged samples (a box car smooth). Write some MATLAB code which smooths the original data – use an index or logical array to exclude the parts of data affected by the noise diode and 1 dB step. Smooth so that you end up with an array smoothed which has a few thousand total points (rather than 127500). Make a note of the code you use to do this below.
7. Instead of using a box car smooth using the mean we can use the median, as that is sometimes more robust to outliers than a mean. Repeat the smoothing exercise above, but use a median window filter instead. Why are the results worse in this case for the median than the mean?
8. Plot the data from the two channels of smoothed data on the one figure. Save this to figure to a JPG or PNG file called exercise5 plt1. Send this plot to your demonstrator at the end of the exercise.
9. If you look at your smoothed data you will see that the off-source part of the scan is not completely flat. We don’t want these baseline variations affecting our Gaussian fitting, so prior to undertaking a non-linear fit we want to remove these. This can be done using the polynomial fitting we developed in the previous exercise. Fit a quadratic baseline to your smoothed data, taking care to exclude the data containing the source, subtract this baseline fit from your smoothed data and create a new array called final. Plot the data from the two channels of smoothed, baseline data on the one figure. Save this to figure to a JPG or PNG file called exercise5 plt2. Send this plot to your demonstrator at the end of the exercise.
10. Repeat this exercise for the data in exercise5b.mat. Plot the data from the two channels of smoothed, baseline data on the one figure. Save this to figure to a JPG or PNG file called exercise5 plt3. Send this plot to your demonstrator at the end of the exercise
2 Non-linear fitting
1. We are now in a position to fit a Gaussian profile to our data. The MATLAB routine we will use to do this is called FM in search. For non-linear fitting it is typically not possible to be certain that the fit returns the global minimum – the parameters space is too large to test all possible values. The concern with non-linear fitting is always that the “best” fit is at a local minimum, rather than the global minimum. There are a variety of non-linear fitting techniques which use different approaches to try and overcome this, but a common factor is that all these methods work better when you pass them a good “initial guess”. Read through the MATLAB help/documentation on min search.
2. Download the fgp.m, gauss profile.m. The code for the fgp (fit gaussian profile) function is relatively straight forward, it is
The first line makes sure that the data to be fitted with a Gaussian pro- file is a column vector. The second line produces an array containing the difference between the data and the passed Gaussian profile and the final line calculates the “error”, a single parameter that attempts to estimate the goodness of fit. This function is passed to FM in search.
3. Have a look at the code for gauss profile. What is the factor parameter in the code?
4. The syntax for calling FM in search is res = FM in search(@fgp,guess,optimset,t, final). Where guess is an array with 3 elements containing an initial guess for the parameters of the Gaussian profile – the amplitude, the FWHM and the centre position. final is the data to fit the Gaussian profile to and t is a “time” array – it can simply be an array of length(final) containing a list of the numbers from 1 to length(final).
5. Use FM in search to fit Gaussian profiles to each channel for each of the two telescope scans you have. Produce plots showing both the final data arrays and the fitted profiles on one figure for one of the channels for one of the scans. Save this figure to a JPG or PNG file called exercise5 plt4. Send this plot to your demonstrator at the end of the exercise.
6. Record the parameters of the fitted Gaussian profiles for each scan/channel below
7. Compare the FWHM of the antenna beam for each of the two scans (which represent scans in two orthogonal directions). How does the measured FWHM compare with theoretical expectations and How close to circular is the antenna beam pattern as judged from these two scans? In order to compare with theory you will need to know the frequency of the observations (8.4 GHz and 2.34 GHz, for the first and second columns respectively) and also the scanning rate of the telescope (to convert from time into angle). These scans were recorded at a scan rate of 1 degree per minute.