Hello Radio: How to Listen to FM Radio Using LimeSDR

LimeSDR is an open source software defined radio which can support a continuous frequency range from 100KHz to 3.8GHz. Many wireless systems such as 2/3/4G, Wifi and FM radio are covered in this frequency range, which makes LimeSDR a fairly useful tool.

To test a software defined radio, usually the first experiment to run is to listen to FM radio. FM radio is the equivalent of Hello World in the universe of the software defined radio. Before starting to describe how to run FM radio test, we have a list of HW/SW needed:

1. LimeSDR board with antenna connected

2. FM transmitter which pairs with a smart phone

3. Lime Suite GUI (version 20.10.0)

4. GNU Radio (version 3.7.11)

The test is done using Ubuntu 18.04.4 LTS.

Step 1: LimeSDR configuration file generation

In the test, we will listen to a FM radio channel set at 89.0MHz. Therefore, LimeSDR needs to generate a tone at 89.0MHz, multiply this tone with the input radio signal, and down-convert the ratio signal to the baseband for further processing. To generate this configuration, in Lime Suite GUI -> SXR, Freq needs to be set as 89.0MHz, then calculate and tune. By pressing the "Save" button at the top left of the GUI, the configuration file will be save as fm_receiver_89Mhz.ini for future usage.

Step 2: GNU Radio set up

GNU Radio is a popular tool for software define radio. To make GNU Radio work together with LimeSDR, one needs to install both GNU Radio and gr-limesdr, a LimeSDR library for GNU Radio. After that, FM_receiver.grc file from /gr-limesdr/examples folder can be loaded into GNU Radio. A few changes need to be made before starting to run it. First, device serial #, .ini file and RF Frequency need to be modified in the LimeSuite Source (RX) module. Device serial # is the serial number for the LimeSDR device and it can be read from Lime Suite GUI; The .ini LimeSDR configuration file has been generated in step 1; RF freq is to be set as 89.0MHz. Baseband Freq in the WX GUI FFT Sink module also needs to be updated to 89.0MHz. Then by pressing the green triangle "Execute the flow graph" button, GNU Radio starts to run. 

Step 3: Run Test

However, you often only hear noise even if the frequency tuned to is a valid local radio channel. The issue is likely due to antenna setting. A rule of thumb is that the antenna size should be in similar dimension as the wavelength of the wireless signal received. For FM radio at 89.0MHz frequency, the wavelength is 300,000,000/89,000,000 = 3.37 meter. But the antenna bought together with LimeSDR is only 4.5 cm long. The antenna size is too small for FM radio. Based on the antenna length, it is designed for multi-GHz freq signal instead of for FM radio which ranges from 88MHz to 108MHz. There are at least two ways to fix it. The first way is to use a different antenna which is much longer. The second way is to move LimeSDR closer, much closer, to the radio transmitter or the other way around. We use the second way. We bought a FM transmitter which can be plugged into the smart phone and tune it into 89.MHz. Then by putting FW transmitter almost right next to the receiver antenna (shown in the second photo below, the first photo shows the my LimeSDR/antenna and their Lego housing), we can now hear the music. 

Task accomplished!

Signal Processing Magic (7) -- Doppler Shift, Rayleigh fading and Jakes model

Doppler shift was discovered by a physicist named Christian Doppler long time ago. Doppler shift is used to characterize the frequency change due to movement. A typical example is that when a fire truck drives toward you, the tone of its siren changes over the time because its speed towards you varies.

Doppler shift for wireless communication is

\[f_d=v/c*f_c\]

whereas v is the speed of the object, c is the speed of light, and \(f_c\) is carrier frequency. For instance, if you sit on a high-speed training moving at 350km/h browsing Internet and your phone is using a 3GHz band, then 

\[f_d = (350*1000/3600)/(3*10e8)*3*10e9 = 972 Hz\]

Based on this equation, moving faster as well as using a band with higher carrier frequency will make Doppler shift larger. Usually the highest Doppler shift is observed in high-speed train. Larger Doppler shift makes channel estimation harder in wireless communication.

Rayleigh fading is a commonly used model for a wireless channel, and it is used when the transmitter can't see the receiver such as in Manhattan island. Rayleigh fading model can be derived from joint Gaussian distribution. Assuming \(x\) and \(y\) are two zero-mean Gaussian R.V. orthogonal to each other, their joint distribution becomes:

\[f_{X,Y}(x,y)=\frac{1}{2\pi\sigma^{2}}e^{-\frac{x^2+y^2}{2\sigma^2}}\]

Rayleigh fading tracks the envelope, \(r\). With \(r^2=x^2+y^2\), we have \(x=r*cos(\theta)\) and \(y=r*sin(\theta)\). To transfer a \(f_{X,Y}(x,y)\) to distribution of \(r\), the determinant of the Jacobian matrix needs to be found:

\[J(r,\theta)=det\begin{bmatrix}\frac{\partial x}{\partial r}&\frac{\partial x}{\partial \theta}\\\frac{\partial y}{\partial r}&\frac{\partial y}{\partial \theta}\end{bmatrix}=det\begin{bmatrix}cos(\theta)&-r*sin(\theta)\\sin(\theta)&r*cos(\theta)\end{bmatrix}=r\]

Therefore,

\[f_{R,\Theta}(r,\theta)=\frac{J(r,\theta)}{2\pi\sigma^{2}}e^{-\frac{r^2}{2\sigma^2}}=\frac{r}{2\pi\sigma^{2}}e^{-\frac{r^2}{2\sigma^2}}\]

After integrating over \(2\pi\) to remove \(\theta\), Rayleigh fading model is:

\[f_{R}(r)=\frac{r}{\sigma^{2}}e^{-\frac{r^2}{2\sigma^2}}\] 

for \(r\geq0\).

In reality, fading channel is always correlated in time. That means the wireless channels your phone sees in this moment and the next moment are always more or less similar. Fourier transmission of wireless channel's time-domain auto-correlation becomes its Power Spectral Density (PSD). For Rayleigh fading with vertical antenna and signal coming from all angles in uniform distribution, it has a famous PSD called Jakes model. Jakes model also depends on Doppler shift, and it is:

\[S(f)=\frac{1}{\pi f_d}\frac{1}{\sqrt{1-(f/f_d)^2}}\]

when \(|f|\leq f_d\) and \(S(f)=0\) otherwise. Jakes model has a famous U shape and the width of the "U" depends on Doppler shift, \(f_d\). When \(f_d=0\), Jakes model is an impulse and it means Rayleigh channel never changes.

Low Light Image Combining Using Python (3) -- Running Even Faster

In this blog, we will discuss how to make lower light image combining run even faster. 

Image alignment is the most time-consuming stage of low light image combining. By using a method similar to what has been proposed in HDR+ project page (https://www.timothybrooks.com/tech/hdr-plus/), we introduce decimation to image alignment and it brings big saving in calculation time.

This is how decimation helps: as discussed earlier, assume block size is B, search range in both horizontal and vertical S, the complexity of image alignment is proportional to B*S*S. For a 4k x 3k image and search range equal to 40 ([-40,40] with each even position), the complexity for matching each pair is 19.2G operations, which is enormous. To save MIPS, we can divide this process to two steps:

Step 1: decimate the image by 16. After decimation, it becomes much easier to align images with the cost of reduced accuracy. The complexity of step 1 now becomes B/16*S*S.

Step 2: we are not done yet. Decimating the original image by 16 means aggregating 16 pixels of the original pixel to one. Now it is time to find out exactly which one of these 16 pixels provides the best alignment so that there is no sacrifice on accuracy. But the search range is drastically reduced. We now align with original image with size B but less range with new S' = 4. S' = 4 is sufficient to find one out of 16 pixel candidates.

Breaking the original image alignment to two steps reduces MIPS quite a bit: instead of B*S*S = 19.2G, now we have B/16*S*S + B*4*4 = 1.39G, a reduction of 93% of complexity. 

By our experiment, the new execution time of the whole image combining process becomes 9.7s without any help from multicore, which is a roughly 1/4 of the baseline. The sample code can be found here: https://github.com/legendzhangn/blog/tree/master/lowlight_image_combine_even_faster

Low Light Image Combining Using Python (2) -- Running Faster

In a previous blog, a Python script to combine several low light raw images for improved SNR has been introduced. But that script runs rather slow. To combine four images together, it takes 155.8s in my PC which uses Intel i5-8400 @ 2.8GHz x 6 core with 16GB memory. How to run it faster? After some attempts, we are able to cut the time to roughly 1/8 of the baseline mainly with two tricks: 1) using decimated image for alignment; 2) processing in multi-core. The new processing times are:

Baseline: 155.8s
After using decimated image for alignment: 40.6s
After using decimated image for alignment + processing in multi-core: 21.1s

Now let me show you how to do it:

Using decimated image for alignment

Image alignment is the most computationally intensive part of the image combining pipeline. In fact, the time spent on image alignment is more than 90% of the total. The operation of image alignment is shown below. During image alignment, the candidate block is swiped though the target block in both horizontal and vertical directions. Assume block size is B, search range in both horizontal and vertical are S, the complexity of image alignment is proportional to B*S*S. For a 4k x 3k image and search range equal to 40 ([-40,40] with each even position), the complexity for matching each pair is 19.2G operations, which is enormous.

One observation is that due to the structure of Bayer pattern, we only calculate candidate of even number of pixel shift such as 0/2/4 etc. Then a natural thought is to decimate both candidate and target by 2. This will bring ~4x acceleration. Instead of B*S*S complexity, it is now B*S*S/4. This explains why processing time is reduced from 155.8s to 40.6s after decimation. Performance wise, this decimation means that instead of using pixels of all colors, only one out of four pixels with green color will be used for alignment. However, the quality of combining seems to hold, which indicate that there are enough pixels left to guarantee the quality of image alignment. The Python code change is below, ":2" is the delta:

Old code:

candidate = rgb_raw_image_candi[f,row_start+row_offset+boundary_adjust[row, col, 0]:row_end+row_offset+boundary_adjust[row, col, 1],
col_start+col_offset+boundary_adjust[row, col, 2]:col_end+col_offset+boundary_adjust[row, col, 3]]
target = rgb_raw_image[row_start+boundary_adjust[row, col, 0]:row_end+boundary_adjust[row, col, 1],
col_start+boundary_adjust[row, col, 2]:col_end+boundary_adjust[row, col, 3]]

New code:

candidate = rgb_raw_image_candi[row_start+row_offset+boundary_adjust[row, col, 0]:row_end+row_offset+boundary_adjust[row, col, 1]:2,
col_start+col_offset+boundary_adjust[row, col, 2]:col_end+col_offset+boundary_adjust[row, col, 3]:2]
target = rgb_raw_image[row_start+boundary_adjust[row, col, 0]:row_end+boundary_adjust[row, col, 1]:2,
col_start+boundary_adjust[row, col, 2]:col_end+boundary_adjust[row, col, 3]:2]

Multi-core processing

Our baseline script use single thread to process. By distributing tasks to multiple threads, we expect the running time to be shorter. Introduction of Python-based parallel processing can be found here. Our task is to align three images with the base image. Therefore, it can be divided into three sub-tasks which are to align each image with the base. In this way, these three sub-tasks can be run independently which brings the best parallel processing gain. Since the sub-tasks are independent, we use the most basic parallel processing module of pooling. What pooling module does is to assign each sub-task to a thread. image_align is the function for sub-task execution and image_input is a list with each element to be input images for a sub-task.

Pooling code:

with Pool(6) as p:
    image_output = p.map(image_align, image_input)

Due to the overhead of parallel processing, the processing time is not cut to 1/3 of single thread. However, parallel processing does reduce the running time by ~20s (40.6s -> 21.1s). In general, when each sub-task is more computationally heavier, there is more time saving by parallel processing.

With both alignment with decimated image and multi-core processing, the final processing time is reduced to 21.1s from 155.8s. My code can be found here 

Low Light Image Combining Using Python (1)

Low light images are photos taken while it is dark. Low light images are usually noisy since there is not enough light coming in. To reduce the noise, a commonly used way is to combine multiple images. For example, such a method has been used in Google pixel camera night mode. Here we will show a Python script which does simple low light image combining.

The pipeline includes three steps: raw image alignment, raw image combining, post processing.

Raw image is the direct output of the image sensor. In normal photo, each pixel has three color channels as R/G/B. But in raw image, each pixel only has one color channel. Bayer pattern is widely used in raw image, and it is also used in this experiment. The color pattern used here is shown below which is one kind of Bayer patterns. As one can see, out of 24 pixels, half of them are green, 25% is red and 25% is blue. The reason that green channel occupies the most pixels is because human eye is the most sensitive to green color.

   0 1 2 3 4 5
0 G R G R G R
1 B G B G B G
2 G R G R G R
3 B G B G B G

What raw image alignment does is to align two raw images. Since multiple images are taken at different times, due to hand movement, there will be small shift of the image. Without alignment, adding images together will create blurry photo. Our method for image alignment is straightforward. To align two images, we will divide both images to blocks with equal dimension. By calling one image "target" and the other image "candidate", we will swipe the candidate block through the target block in both x and y axis search directions. The scope of the swipe is called search range. For each swipe position, the L1 distance will be calculated which is the delta of the pixel values of the candidate and target blocks. Out of all swipe positions, the one with smallest L1 distance will be selected. For more details of this operation, you can refer to the Python source code here.

The reason that image alignment should be done at the very beginning of the image pipeline on raw image is for both performance and computation time. In term of performance, additional processing of the signal at the later stages such as interpolation/correlation may degrade the performance of alignment. In term of computation time, raw image has one color in a pixel but at later stages after interpolation, there will be three colors in each pixel. Calculating alignment on one color channel is faster than calculating that on three channels. Due to the structure of Bayer pattern, we only calculate candidate of even number of pixel shift such as 0/2/4 etc.

The step of raw image combining is to add candidate and target blocks together. Since they are already aligned, adding them together will enhance the signal while reducing the noise. Here we apply equal weight to each image. However, performance could be improved by applying different weights.

Between raw image and the photo we usually see, post processing needs to be done. Libraw library is used for post processing. The major steps of post processing is shown below:

Next we will show the outcome of combining. Raw images are obtained using Pro camera mode of a Samsung Galaxy S9 device. The camera settings are ISO = 800, 1/4 shutter speed, and F1.5. S9 has two aperture modes: F2.4 and F.15. F1.5 is a larger aperture for low light photo. Figure below shows the difference before (left) and after (right) the combining. The image quality is clearly improved. To get the right figure, we combine eight low light images. Bilateral filtering can be used to further enhance the image quality.

Source code of the Python script and raw image files can be found here: https://github.com/legendzhangn/blog/tree/master/lowlight_image_combine

LibRaw Post Processing Pipeline Translated Into Python (2) -- Low Light Photo

In a previous post, we introduced a Python script which emulates the pipeline of LibRaw. LibRaw is a commonly used library for raw image conversion.

Then we found that compared with LibRaw, our Python script does not work well for raw image taken under low light. While LibRaw still shows the objects albeit quite noisy, our original Python script does not show anything. After some digging, it is found that LibRaw has an auto bright function which is not available in the previous Python script. As shown below, for low light image, the difference made by using auto bright function is quite stark.




Thereafter, our Python script gets updated and now it supports auto bright as well. The main difference between low light image and normal light image is how Gamma mapping is done. For low light image, Gamma curve peaks much faster than for normal light image. The exact Gamma curve of low light image depends on its histogram. The updated script can be found here, and it shows how this new Gamma curve is generated.

An Android App for Photo Brightening

I published an Android app used for photo brightening named Brighter Photo. Using a low light photo as the input, the app can post process it and make it look brighter. It is used only for post processing but not as a camera. The app is available at Google Play store and the link is here

Figure below shows what the app can do. The left half is a photo taken at night and the right half is the same photo after processed by Brighter Photo app.

To use the app, the GUI has three buttons: Load, Save and Setting. "Load" button can load photo from Android gallery. After loading, post processing immediately starts. Post processing takes 10-30 seconds depending on the image size. It runs slower than some of the similar apps but produces better results. "Save" button allows to save the results. Note that to enable saving, the storage permission needs to be enabled in Android setting for this app.

This Youtube video shares how to use the app.

Hope that this introduction is useful and you can enjoy using this app. Any feedback is welcomed.

A Simple Example of Image Processing Using Java

In this post, we share a simple example of image processing using Java. This example benefits from Java tutorials available from Internet.

This example has three steps: 1) read the input image; 2) image processing; 3) display the processed image.

For step 1, an image is read into image buffer. 

    File file = new File("Lenna.png");
    BufferedImage image = null;

  // Read image
    try
    {
        image = ImageIO.read(file);
    }
    catch (IOException e)
    {
        e.printStackTrace();
    }
    System.out.println("done");

Step 2 has a simple operation of image processing. It does 2D convolution filtering for the image using filt3x3 function in the ConvolutionMatrix class.

    // 2D convolution for the image
    double[][] config = {{1,2,1}, {0,0,0}, {-1,-2,-1}};
    ConvolutionMatrix imageConv = new ConvolutionMatrix(3);
    imageConv.applyConfig(config);
    image2 = imageConv.filt3x3(image, imageConv);

Step 3 displays the modified image

    // Display the modified image
    ImageIcon icon=new ImageIcon(image);
    JFrame frame=new JFrame();
    frame.setLayout(new FlowLayout());
    frame.setSize(image.getWidth(),image.getHeight()); //Window.setSize(int width, int height)
    JLabel lbl=new JLabel();
    lbl.setIcon(icon);
    frame.add(lbl);
    frame.setVisible(true);
    frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);

Java-based image processing is used in many places. Hopefully some people can benefit from this tutorial. The code can be found from here

LibRaw Post Processing Pipeline Translated Into Python

LibRaw is a widely used open source library for raw image conversion. In many people's opinion, the most valuable part of LibRaw library is that it can handle various types of raw image formats, which is not commonly available to other libraries. But people are also interested in its post processing pipeline, which converts a interleaved RGB dot map to a colorful picture by using procedures such as white balance adjustment, demosaicing, gamma mapping, etc. People want to understand LibRaw's post processing pipeline but the code can be hard to deciphered. That is why I wrote a Python script which appears to match with LibRaw's post processing pipeline quite well. The hope is that by reading this Python script, it is easier for folks to understand what actually happens in LibRaw.

I should admit that some corners have been cut in this Python script. For example, the output of the script is only half the size of the input, which means no need for me to write a full-fledged demosaicing which matches the size of the input. However, what I found over the time is that quite often, the most puzzling part of post processing to people include me is how inputs from different color channels are scaled and mixed together. The scaling and mixing parts are covered in this Python script.

Here is the approach we take when writing this script: the input of this script is a raw .dng image. The image is taken using Samsung Galaxy S9 phone. The decoding before post processing is done using Python rawpy library. Python rawpy library is LibRaw wrapped in Python. After decoding, post processing is done using two parallel paths. The first path is to use postprocess() function provided by rawpy. Under the hood, it calls dcraw_process() from dcraw library written by Dave Coffin. The first path is used for reference. The second path is my own pipeline written in Python. The second path takes rawpy decoding outputs as inputs including raw image and some parameters. At the end of the script, we compare the reference and reconstructed image by my code and find them to be almost the same.

Now let me introduce the post processing code I wrote. It includes three main steps: 1) scale the input by white balance, 2) mix R/G/B by using color matrix, 3) Gamma scaling

In step 1, raw image containing RGB info is scaled by white balance (wb_scale). Since Bayer pattern is used for raw image input, there are two green channels (color4[:,:,1] and color4[:,:,3]). The white balance scale comes from the parameter of camera_whitebalance in rawpy decoding output. One can use another parameter of daylight_whitebalance but it makes the final output image yellowish. Since the next step of pipeline only needs three inputs, there is minor step called mix green which mixes two green channels as: color4[:,:,1] = (color4[:,:,1] + color4[:,:,3])/2.

color4[:,:,0] = raw_py.raw_image[0::2,1::2]*wb_scale[0]
color4[:,:,1] = raw_py.raw_image[0::2,0::2]*wb_scale[1]
color4[:,:,2] = raw_py.raw_image[1::2,0::2]*wb_scale[2]
color4[:,:,3] = raw_py.raw_image[1::2,1::2]*wb_scale[3]

In step 2, color_matrix is applied. I believe the purpose is convert camera sensor input to standard RGB image. color_matrix is also among rawpy decoding outputs. While white balance changes from photo to photo, color_matrix is camera's inherent property and it stays the same with the same camera setting.

pic_temp[:,:,0] = color_matrix[0,0]*color4[:,:,0] + color_matrix[0,1] * color4[:,:,1] + color_matrix[0,2] * color4[:,:,2]
pic_temp[:,:,1] = color_matrix[1,0]*color4[:,:,0] + color_matrix[1,1] * color4[:,:,1] + color_matrix[1,2] * color4[:,:,2]
pic_temp[:,:,2] = color_matrix[2,0]*color4[:,:,0] + color_matrix[2,1] * color4[:,:,1] + color_matrix[2,2] * color4[:,:,2]

Step 3 is Gamma mapping. Gamma mapping is a legacy from CRT era to compensate the nonlinear distortion of CRT monitor. But it outlives CRT. The parameter used in Gamma curve is Gamma(0.45, 4.5).

pic_temp = gamma_curve[pic_temp]

Finally we compare the results between reference image by rawpy's postprocess() function and reconstructed image by our code. It turns out that max delta between them is 2 output 255. Thus we claim that these two images are almost the same. The code/image can be found here.

Image White Balancing with Python

Below is a photo I took in the viewing deck of Taipei 101. Taipei 101 has a glass wall which is green color. Thus, all photos taken inside the building through the wall has a green cast including this one.

To remove the cast of green color, we use the technique of image white balancing. First, we need a reference object with known color. Fortunately this photo has clouds (pointed by the arrow). Clouds are good references since they are white. By averaging the pixels of cloud (mean(img_array[400,500:550,:])), the mean of R/G/B is [201.5, 254.9, 253.9]. Since white color has R/G/B channels roughly equal, to make it white, 52 needs to be added to R channel of the whole image. Like DC cancellation in communication, white balancing removes the bias of the signal.

Code below shows how to do white balancing

filename = 'IMG_9254.JPG'
img = image.load_img(filename)
img_array = np.array(image.img_to_array(img),dtype=np.uint8)
img_shape=img_array.shape
plt.figure()
plt.imshow(img_array)

img_array2 = img_array
img_array2 = np.array(img_array2, dtype=np.uint16)
img_array2[:,:,0] = img_array2[:,:,0] + 52 # white balancing
img_array2 = np.clip(img_array2, 0, 255)

The newly generated image after white balancing, img_array2,  is shown below. The color does look more natural. In case that you want to repeat this, code/image can be found here.

Note on Implementation of “Fast Noise Variance Estimation” by J. Immerkær

Noise variance estimation is a fundamental task of image processing. Among various approaches of noise variance estimate, fast noise variance estimation method proposed by J. Immerkær in [1] stands out because of its good balance between complexity and accuracy.

Based on J. Immerkær's method, to estimate the standard deviation of image noise \(\sigma\), we can use equations below:

\(F=
\begin{bmatrix}
1 & -2 & 1\\
-2 & 4 & -2\\
1 & -2 & 1\\
\end{bmatrix}
\)

\(\sigma=\frac{\sqrt{\pi/2}}{6(W-2)(H-2)}\sum_{I}|I(x,y)*F|\)

\(F\) is the high-pass filter. To estimate noise variance, filter \(F\) will be applied on the whole image \(I(x,y)\) with convolution. The sum of absolute value of convolution results is normalized by width of height of the image, which is \((W-2)(H-2)\). "-2" counts for the boundary effect. The summation is also scaled by \(\frac{\sqrt{\pi/2}}{6}\), which is explained next.

Assume a 3x3 block of \(I(x,y)\) with each element of \(x_{i}\) to be a Gaussian random variable with mean \(\mu\) and standard deviation \(\sigma\)
\(I=
\begin{bmatrix}
x_{1} & x_{2} & x_{3}\\
x_{4} & x_{5} & x_{6}\\
x_{7} & x_{8} & x_{9}\\
\end{bmatrix}
\)
Then \(|I*F|=|x_{1}-2*x_{2}+x_{3}-2*x_{4}+4*x_{5}-2*x_{6}+x_{7}-2*x_{8}+x_{9}|\)
Since \(E(|I*F|^2)=36\sigma^2\), \(|I*F|=6|s|\) where s is Gaussian random variable with mean 0 and standard deviation \(\sigma\). We can derive the distribution of \(y=|s|\) and it is \(\frac{2}{\sqrt{2\pi}\sigma}e^{-\frac{y^2}{2\sigma^2}}\)  for \(y>=0\) and 0 for \(y < 0\). <0 .="" br=""> 

Since \(E(|s|)=\sqrt{\frac{2}{\pi}}\sigma\), it explains why the summation is also scaled by \(\frac{\sqrt{\pi/2}}{6}\).


[1] J. Immerkær, “Fast Noise Variance Estimation”, Computer Vision and Image Understanding, Vol. 64, No. 2, pp. 300-302, Sep. 1996

Google Cloud Service GPU Profiling: K80/T4 vs GTX1060

Recently I tried training neural network using Google cloud service. The user experience is in general good. The 300$ free credit also helps. To support deep learning task, one can create virtual machines in Google cloud service which contain both CPU and GPU. There are a few options of GPU ranging from K80 to V100 depending on how heavy the computation task is and user's affordability. I compared the performance of two GPUs: K80 and T4 with my local machine. Here are the results:

The benchmark used is a Keras example based on MNIST data set. The metric for profiling is the duration used for one epoch of training. The table below shows the profiling results:

GPU	CPU	Time
GTX1060 6GB GDDR5	Intel 6-Core i5, 16G RAM	35s
Tesla K80, 12GB GDDR5	4vCPU, 26G RAM	8s
Tesla T4, 16GB GDD6	4vCPU, 26G RAM	4s

GTX1060 is the GPU used in my local machine. The profiling results show that for this MNIST benchmark, the time used by K80 is about one fourth and the time used by T4 is about one ninth of that of my local machine. Note that when CUDA/Tensorflow libraries are not set up correctly, the computation might be done in CPU and the computation time will be increased drastically. For example, the duration of training one epoch will increase from 8s to 84s in case the computation is done in CPU not GPU. Thus when computation time is unexpectedly long in Google cloud service, please confirm that the computation is performed indeed in GPU.

Signal Processing Magic (6) -- Half Band Filter

Half-band filter is a FIR filter often used in decimation by 2 filtering operation. Benefit of this filter is that half of the filter coefficients are 0's, which reduces the complexity of filtering operation. The name of "half-band" comes from that the filter's cut-off frequency is roughly fs/4, thus the filter's pass band is around half of the total.

To design a half-band filter, it starts from a sinc function as \[s[n]=\frac{sin(\pi*n/2)}{\pi*n/2},n=-N,-(N-1),...,N\]
We know that when N approaches infinite, freq response of this sinc function is rectangle ranging from -fs/4 to fs/4. But in real world, N is always a limited number, and when that occurs, the filter does not look great. Ripples occur in pass band and rejection in stop band is often not good enough either.

How to improve this? Windowing method comes for rescue. As introduced before, windowing method enhances the pass band (main lobe) and suppresses the stop band (side lobes). Assuming window is w[n], after adding window, the new filter becomes
\[h[n]=s[n]*w[n]\]
Since half of \(s[n]\)'s coefficients are 0's, half of \(h[n]\)'s coefficients are 0's too.

Matlab script to generate a 21-tap half-band filter with Hamming window is shown below:

N = 10;
x = [-10 : 10]/2;
s = sinc(x)
w = hamming(2*N+1);
h = s.*w.';

figure;
subplot(2,1,1);stem(s);title('Sinc');
subplot(2,1,2);stem(h);title('Half-band');

[H, F] = freqz(s, 1, 'whole', 2^15);
[Hf, Ff] = freqz(h, 1, 'whole', 2^15);

figure; hold on;
plot(F(1:end/2)/pi/2, 20*log10(abs(H(1:end/2))));
plot(Ff(1:end/2)/pi/2, 20*log10(abs(Hf(1:end/2))));
grid;
ylabel('Magnitude (dB)'); xlabel('Freq (fs)');
legend('Sinc','Half-band');

Filter coefficients for Sinc and half-band filter

Freq response for Sinc and half-band filter

TensorFlow Android App Debugging (3) -- Keras to TensorFlow Model Conversion

Deep learning model can be programmed using different libraries. Among them, Keras is one of the most easy to use. However, network model generated using Keras can''t be directly applied to Android app. First it needs to be converted to TensorFlow. A conversion tool is provided here to automatically convert Keras model to TensorFlow. The instruction of this tool is straightforward. It converts .h5 (Keras model) to .pb (TensorFlow model).

    python keras_to_tensorflow.py 

        --input_model="path/to/keras/model.h5" 

        --output_model="path/to/save/model.pb"

For verification purpose, we can use tool named Tensorboard to visualize the generated model. As the first step, the model will be read into a folder by script read_pb_model.py. The default folder name is tf_summary.

python read_pb_model.py model.pb

Then tensorboard will convert the folder's content to a graph, and display it in a webpage such as http://nan-System-Product-Name:6006. Below is an instruction of tensorboard 

tensorboard --logdir ./tf_summary

An example of generated graph is shown below. The graph displays the network topology. It also shows the name of the node. It is important to know the names. In Android app, when one needs to access certain part of the network model, he has to use the name as reference.

For CV2 Color to Grayscale Image Conversion

The task is simple: we want to verify the formula used in CV2 (OpenCV 2.0) for converting color image to grayscale. The formula to convert from RGB to gray can be found here: Y = 0.299*R + 0.587*G + 0.114*B

We use the following code to verify the formula. First we read in a color image as both color (img_color) and grayscale (img_gray). Then we use the equation above to convert the color image to gray (img_color2gray). As the final step, we compare img_gray and img_color2gray. This comparison allows to verify the equation. The output of the equation is floating number. We use rounding (np.round) to convert float to integer. np.round(1.499) = 1 and np.round(1.5) = 2. We have tried other quantization methods such as floor or ceil but they don't behave as well.

img_name = 'lena_color_512.tif';
img_color = cv2.imread(img_name, cv2.IMREAD_COLOR)
img_gray = cv2.imread(img_name, cv2.IMREAD_GRAYSCALE)

shape = img_gray.shape
img_color2gray = np.zeros(shape, dtype=np.int32)
img_color2gray[:,:] = np.round(0.299*img_color[:,:,2]+0.587*img_color[:,:,1]+0.114*img_color[:,:,0])
img_color2gray = np.clip(img_color2gray, 0, 255)

The comparison result is this: img_gray and img_color2gray are very similar but not identical. Out of 512*512=262144 pixels in total, 352 pixels don't match (%0.13 of total). The dots in the figure below are the unmatched pixels. This leaves a question of how exactly this conversion is done in OpenCV SW. My source code can be found here.

Signal Processing Magic (5) -- Windowing

In discrete-time signal processing, if received signal is a tone, its spectrum looks like a sinc function (see this). The reason is by cutting digital signal to finite length, it is like imposing a rectangular window in time domain. And rectangular window in time domain corresponds to sinc function in freq domain.

But some engineers thought side lobes in sinc function are too big. In order to make them smaller, they invented windowing method, which is to impose a window on the time-domain signal and it can suppress side lobes in the frequency domain. Since then people have invented different types of windows (Bartlett, Hanning, Hamming, Blackman etc. see Sec 7.2 in Discrete-time signal processing written by Oppenheim&Schafer). Hamming and Hanning windows are shown below. Assuming the digital signal has M samples, to apply windowing, one should generate a window lasting M samples and then multiple window with digital sample by sample. Each window is well defined mathematically. For example, Hamming window for M samples is defined as:
\[w[n]=0.54-0.46cos(2\pi n/M)\]

Next we compare two signals, one with windowing and one without. The signal is a 100KHz tone with sampling rate of 1MHz. The plot below shows that by adding Hamming window, the side lobe is suppressed by around 20dB but main lobe becomes wider. Windowing function works like a magical rolling pin. It drives the signal energy from side lobes to the main lobe.

In the next plot, instead of one tone, we show the spectrum of two tones, one at 100KHz and another at 300KHz. After applying windowing, these two tones become more distinctive.

Python Implementation for BM3D Denoising of Color Image

BM3D denoising is a popular method used for removing image noise. The original BM3D paper [1] proposed the algorithm for gray scale image. Then in [2], the method is extended for color images.

Python script for grayscale image BM3D can be found here by liuhuang31. But Python script for color image BM3D can't be identified after searching in Internet. Note that this Python implementation supports BM3D for color images but the kernel part is written in C by Marc Lebrun. By modifying liuhuang31's script, we provide a Python implementation for color BM3D here. In this implementation, all codes including algorithm kernel part is written in Python.

Below is an example including noisy image and images after step 1 and step 2 of denoising processing of BM3D. Noisy image is generated by adding random noise to reference noise-free image. PSNR is peak SNR calculated between each individual image and the reference noise-free image.

Noisy image (PSNR = 22.41dB)

After step 1 of removing noise (PSNR = 29.36dB)

After step 2 of removing noise (PSNR = 30.13dB)

While our Python script is an extension of liuhuang31's work, we should note that we have made a few corrections:

Correction 1:
In the original code, delta of two images is calculated as below. But since img1 and img2 are both uint8 type, the calculation result is wrong.

D = numpy.array(img1 - img2, dtype=numpy.int64)

In the modified code, img1 and img2 are first typecast to int64 and then delta between them is found.

D = numpy.array(img2, dtype=numpy.int64) - numpy.array(img1, dtype=numpy.int64)

Correction 2:
In the original code, to find boundary, it uses the following code. But last line has a typo. shape[0] should be shape[1].

    if LX < 0:   LX = 0
    elif RX > _noisyImg.shape[0]:   LX = _noisyImg.shape[0]-_WindowSize
    if LY < 0:   LY = 0
    elif RY > _noisyImg.shape[0]:   LY = _noisyImg.shape[0]-_WindowSize

The modified code is:

    if LX < 0:   LX = 0
    elif RX > _noisyImg.shape[0]:   LX = _noisyImg.shape[0]-_WindowSize
    if LY < 0:   LY = 0
    elif RY > _noisyImg.shape[1]:   LY = _noisyImg.shape[1]-_WindowSize

Correction 3:
The last one is probably the trickiest one. To calculate Wiener filter in step 2, noise variance is sigma^2. To match with noise variance, signal power should be normalized by the count of similar blocks. This normalization is not in liuhuang31's original code. Other BM3D code such as Marc Lebrun's code also uses this normalization.

                Norm_2 = numpy.float64(tem_Vct_Trans.T * tem_Vct_Trans)
                m_weight = Norm_2/Count/(Norm_2/Count + sigma_color[ch]**2)

[1] K. Dabov et. al., Image denoising by sparse 3D transform-domain collaborative filtering, IEEE Trans. Image Proc., Vol 16-8, Aug 2007
[2] K. Dabov et. al, Color image denoising via sparse 3D collaborative filtering with grouping constraint in luminance-chrominance space, ICIP 2007

TensorFlow Android App Debugging (2) -- Dump Out Image/Text Files

During Android app debugging, dumping out data from Android to PC to crucial for the following analysis process. In the previous post, we said that Android logging function is a useful tool. But sometime logging only is not sufficient. For example, to save information of an image by logging is time-consuming, prone to packet loss and unproductive. A better way is to save the image into internal memory. In this post, we provided some examples of how to save images and text files in Android.

To save images (modified from the codes found here). bitmap is saved into .JPEG image files with lossless compression.

    // Save bitmap
    String extStorageDirectory = Environment.getExternalStorageDirectory().toString();
    OutputStream outStream = null;
    File file = new File(extStorageDirectory, "/DCIM/bitmap"+fileNo+".JPEG");
    try {
      outStream = new FileOutputStream(file);
      bitmap.compress(Bitmap.CompressFormat.JPEG, 100, outStream);
      outStream.flush();
      outStream.close();
    } catch(Exception e) {

    }

To save text files (modified from function found here)

  private void writeToFile(String content, final String filename) {
    try {
      File file = new File(Environment.getExternalStorageDirectory() + "/DCIM/" + filename);
      //File file = new File("/DCIM/test.txt");

      Log.i(TAG, Environment.getExternalStorageDirectory() + "/DCIM/" + filename);
      if (!file.exists()) {
        file.createNewFile();
      }
      FileWriter writer = new FileWriter(file);
      writer.append(content);
      writer.flush();
      writer.close();
    } catch (IOException e) {
    }
  }