Why does the frequency in the DFT have to be an integer? I don't know why the result of DFT and FFT in MATLAB like below images..

fs=128;
t=0:1/fs:1-1/fs;
x=cos(2*pi*3.5*t);
X=fft(x);
N=length(x);
n=0:N-1;
f=fs*n/N;
plot(f,abs(X)/N);

If I set the frequency to an integer, the result of above is like this

fs=128;
t=0:1/fs:1-1/fs;
x=cos(2*pi*3*t);
X=fft(x);
N=length(x);
n=0:N-1;
f=fs*n/N;
plot(f,abs(X)/N);

Please teach me why the difference is occured..

There are a few different ways to interpret the math of the DFT. The one I find most suitable to explain this and other properties of the DFT assumes that the continuous Fourier domain was sampled, which causes the time domain to become periodic (just like the signal was sampled, causing the Fourier domain to become periodic).

If your signal contains a sine wave, it has a single frequency. But if a non-integer number of periods of this sine wave is represented by the samples in your discrete signal, then there will be a discontinuity when replicating these samples. This discontinuity contains all frequencies.

The "frequency" parameter used in the DFT, often represented as k, is in relation to the length of the signal. If your signal shows exactly a single period of the sine wave, this sine wave has a frequency k=1. If you see exactly two periods, it has a frequency k=2. If you see one and a half periods, then k would be 1.5, which does not exist.

Draw one and a half period of a sine wave. Now copy that signal multiple times end to end:

x = cos(linspace(0,3*pi,1000));
plot(repmat(x,1,7))
set(gca,'xticklabel',0:7)

You now have a signal where the lowest frequency corresponds to the repetition of the original signal, the period length is the length of the original signal. This is k=1. The next frequency is twice this lowest frequency (period is half). You will notice if you study the graph above that there is nothing that repeats with a frequency of 1.5 times the lowest frequency (period is 1.5 times the original signal). Such a frequency does not exist because of how the signal was constructed.


is $F(x, n)=x(n)u(n)$ time invariant? My book with answers to proakis and manolakis book says that $F(x,n)=x(n)u(n)$ is time invariant while my calculation shows that it is time variant. I have repeated them many times and I am mostly convinced that book is wrong. Can somebody please check my calculations?

let $u(n)$ be a Heavside step function. First we shift output signal:

$$F(x, n-k)=x(n-k)u(n-k)$$

Let $x'(n)=x(n-k)$ be a shifted input signal. Lets calculate response to shifted input signal:

$$F(x', n)=x'(n)u(n)=x(n-k)u(n)$$

since $F(x, n - k) \neq F(x', n)$ this system is time variant. Is my method correct?

The answer to this (and your other similar question) is most likely found in the somewhat unusual notation used for the system.

In ordinary terms, if you just take your system as a map between input and output signal, then you're absolutely right, it would not be time invariant.

However, Proakis defines the system as a map from an input signal and a time coordinate to the output signal. If you want, the time coordinate itself can be understood as a second signal that just happens to contain the time index at each time index. Consequently, the time shift operator would act on this system definition by both shifting the proper signal and the time coordinate, so that the the system turns out to be time invariant.

So let me summarise: The system map $F(x,n)=x(n)u(n)$ is time invariant whereas the map $F(x)=x(n)u(n)$ is not time invariant.


Scrambler generator polynomial with initial state A length-$127$ frame scramble uses the following generator polynomial

$$G(D)= D^7+D^4+1$$

If the all-ones initial state is used, how can one obtain the $127$ bit sequence? I am not interested in the sequence itself or computer algorithm but an understanding of how to obtain in.
My understanding is that $G(D)$ is the $\mathcal Z$-transform of the impulse response $\textrm{1001001}$. Is this convolved with the initial state $\textrm{1111111}$?

This generating (check if it's a primitive polynomial - those are usually used) polynomial representing a simple Linear Feedback Shift Register (LFSR). The circuit is given below, the gate is a XOR gate (an addition over GF(2) ). You can follow the output of this shift register very easily


How to find period/envelope that is larger than sample size? Mercury's distance from the solar system's barycenter (in kilometers)
can be approximated as:

5.9115960588705115e7 + 1.1573608483842954e7*Cos(2.086944367456105 - 0.0029751428085695573*t)

based on 8 years worth of hourly "samples" from HORIZONS. (t is in
hours; I forget what value corresponds to t=0, but it should be
irrelevant to my question).

The residuals after this approximation look like this:

I could add additional Fourier coefficients, but there appears to be a
period or envelope that is larger than the sample size.

Standard Fourier analysis will never find this period, and continuous
Fourier analysis is fairly inefficient.

How can I find this period/envelope?

Perhaps more to the point: what's the simplest function that does a
"fairly good" job of approximating this data?

The graph shows a pattern, so there "must be" such a function?

If you can guess the number of Fourier coefficients, then a decomposition method such as MUSIC or ESPIRIT might provide a parametric estimate or eigenvector decomposition, given that guess.


[2D fourier transform]: Most people can't explain this I am confused with 2d magnitude plot of frequency spectra. So we have 2 images, the first one is shown at the top, and the dilated or enlarged version of the white box is shown at the second row.

For the dilated version, in the spatial domain when travelling in the $x$-direction, you see a fast change in intensity, and when you travel in the $y$ direction, the change is slower, so that the $v$ direction frequency component is smaller. In the frequency domain, $x$ maps to $u$, and $y$ maps to $v$.

I do understand that when you enlarge in one domain, you must do the opposite to the other domain.
The thing I don't understand is that in the dilated version, when you travel in $x$, and you still see the same rate of change as in the original image. Since $u$ maps to $x$ in the frequency domain, we expect that the frequency spectrum looks the same in the $u$ direction. However, the result is that, in the $u$ (or the horizontal) direction of frequency domain, its frequency component shrinks; the $v$ direction remains intact.  But we just change its $y$ direction!

Update:
To be more precise, please refer to the area of the red box. In the second figure, it gets shrink.

The scaling property of the Fourier transform tells you that if  $\mathcal{F}(x(t)) = X(f)$, then  $\mathcal{F}(x(at)) = \frac{1}{|a|}X(\frac{f}{a})$, with $a\neq 0$. So when you dilate in time, you contract in frequency, and vice-versa.

Now, for continuous 1D signals, the Fourier transform of a uniform box will be a cardinal sine. If you dilate the box, the cardinal sine will be contracted. This could be coarsely observed from the location of its zeros.

Here, you deal with a 2D discrete signal. So the result should be a discrete (or periodic) version of a cardinal sine. Its spread in one direction can be visually estimated by the horizontal and vertical black lines, since your 2D discrete Fourier transform is separable.

Now look at one horizontal line in your red boxes. The blackest dots that separate the horizontal gray/white bumps have the same location on both images. Now look at several lines. The black dots align vertically. So you see that the vertical lines are spaced in the same manner on the top and the bottom, which is exactly what you expect.

The result in the vertical direction is the controverse. Horizontal black lines are closer to each other in the bottom row.

What is sometimes misleading in Fourier plots is that some features in the frequency domain look orthogonal to those in the spatial domain. You can see more of this effect in An Intuitive Explanation of Fourier Theory.

If you want to gain additional conviction, crop a line from each spatial image, and an horizontal line from the spectrum, and superimpose the plots.


Non radix-2 FFT hardware implementation : zero-padding vs custom algorithms I am an FPGA engineer with some DSP experience, facing the situation of implementing a 1023 points FFT in a spread spectrum demodulator.

As a straightforward solution, I would zero pad to 1024 and use an IP. But I have been given the task to optimize the design : the key parameter here is area, since the sample rate is really low (< 30 Msps) and there is possibility of reusing logic, but latency optimization would be also interesting.

My question is: can anyone refer me to bibliography hinting non-radix 2 FFT implementations comparing performance/implementation constraints tradeoffs?

thanks

The FFT is just a fast method to implement the Discrete Fourier Transform (DFT). It's based on the idea that you can break a large DFT down into a bunch of smaller ones and the combine the results for the final DFT. That requires the DFT length to be broken down into it's prime factors. The more and the smaller the prime factors are, the more efficient the algorithm. That's why a power of 2 is the best choice: the prime factors are small (they are all 2) and there is a lot of them.

1023 is awkward since the prime factors are 3, 11 & 31, so the efficiency gain you can get from an FFT over a straight DFT algorithm is limited. Stating the obvious: a 1024 FFT would be much faster, take a lot less code and data memory and will essentially have the same latency as a 1023 FFT. You better have a VERY good reason to not zero-pad one sample.

For the basics on prime factor FFT see https://en.wikipedia.org/wiki/Prime-factor_FFT_algorithm.


Is it possible to convert from 2D to 3D form? I have captured car positions at different frames from the video.
Suppose car's(e.g.left side car which is coming towards us in the images) centroid  at video frame1 is P(x1,y1) and Q(x2,y2) at video frame4.

Is it possible to represent P and Q points into 3D? so that i can calculate correct pixel distance d(PQ) and finally compute actual distance

Note:u can assume that camera is stationary, placed at height 10 m from ground level .u can also assume any suitable data if u want
also u can refer the next figure also.

[EDITED] Here's how it's done.

Steps:
1. Isolate the Road Divider Part.

Then, using Houghlines, find out the longest lines in Image. Find out the extrema points that cross image boundary. You got the Quadilateral points. I skipped this part by Manually Choosing them. In my case, the width of road at top of image is 10, and at bottom is 60.

Now, for the Destination, we have to Image so that Divider Quadilateral has uniform width of 60, so we save in another array new points, which have same bottom points, but top points are modified as they form a rectangle with height same as image (or bigger, you can adjust as it suits you after running it once and seeing the results), but width 60. Remember, we are transforming our image such that the divider fits into this Rectangular Area. Image will change automatically.

Now, using
Mat TransformMat = getPerspectiveTransform(ipPts, opPts);

of opencv on initial and destination quadilateral Points, we get a Transform Matrix, Which we apply to our Test Image.

warpPerspective(ipImg, opImg, TransformMat, ipImg.size());

You'll see a lot of Area is Cropped to fit the Output in InputImage Size. However you'll notice that Image is Transformed to get parallel Divider, as we Expected. By performing some Translation and Stuff, and taking a bigger Size Output Image, You'll get this. We can never get a perfectly wrapped Output Image, since its size can be very large. But we'll reach some constraint that'll fulfill our task. Here's Last image :


How to Prove a 2D Filter Is Separable? I want to prove that 2D Gaussian filter is separable and we can separate it into two dimensions, my problem is about the size of filters. we should prove that $G(x,y)*I$(where $G(x,y)=$$\begin{bmatrix}0.01 & 0.1 & 0.01 \\0.1 & 1 & 0.1 \\ 0.01 & 0.1 & 0.01\end{bmatrix}$,I is image and $*$ is convolution operator) is equal to $G(x)*I*G(y)$ where $G(x)=$$\begin{bmatrix}0.1 & 1 & 0.1  \end{bmatrix}$ $,G(y)=$$\begin{bmatrix}0.1\\1\\ 0.1\end{bmatrix}$ and I is image. in other words we should prove that $G(x,y)=G(x)*G(y)$ but I don't know how to convolve these filters with different sizes.

Let's have a different perspective on that.

Let's say our 2D Linear Operator is given by the Matrix $ G \in {\mathbb{R}}^{n \times n} $.
Using the SVD Decomposition the operator can be written as:

$$ G = \sum_{i = 1}^{n} {\sigma}_{i} {u}_{i} {v}_{i}^{T} $$

Separable Linear 2D Operator is defined as operator which can be composed by Outer Product of 2 vectors.
Looking at the SVD Decomposition of $ G $ we can conclude that $ G $ is separable operator if and only if $ \forall i > 1 \; {\sigma}_{i} = 0 $ and it is given by:

$$ G = {\sigma}_{1} {u}_{1} {v}_{1}^{T} $$

Usually LPF 2D Linear Operators, such as the Gaussian Filter, in the Image Processing world are normalized to have sum of 1 (Keep DC) which suggests $ {\sigma}_{1} = 1 $ moreover, they are also symmetric and hence $ {u}_{1} = {v}_{1} $ (If you want, in those cases, it means you can use the Eigen Value Decomposition instead of the SVD).

So basically, to prove that a Linear 2D Operator is Separable you must show that it has only 1 non vanishing singular value.

Note: I was not rigorous in the claims moving form general SVD to the Eigen Decomposition yet the intuition holds for most 2D LPF operators in the Image Processing world.


Is HSV Color Space Sufficient for Rudimentary Color ID and Edge Detection I want to, as part of a C++ school project (no OpenCV), work with an RGB image taken from a camera filming smarties (round colored pellets) from above and extract the number of candies in each color.

Is converting the (reasonably lit and clear) image to HSV sufficient to perform edge detection (with Sobel on the V channel) and color identification (with H and S), following noise removal and possibly histogram (or other kinds of) equalization? Or should I stick with RGB color space and do more thorough filtering in order to properly isolate the same color of differently lit smarties.

Hope this makes sense, thank you for your time

The approach seems reasonable.
Indeed doing edge detection in weighted RGB channel is the classic approach (Though you could also employ more advance methods, See Edge Detection on a Color Image).

I think you could achieve great results if you also look specifically for oval shapes then you reduce the chances for false positives.

Color identification in the H and S channels will give most information (Well, up to differences in luminosity of the color).


How can I trace a path around an image consisting of only black and white lines? I'm looking for a way to trace a closed curve (i.e. a long list of pixel values all connected by lines in between) that approximates a bitmap image of black and white lines.

The image starts out as an arbitrary bitmap image. I've been using this one for testing.

Then it's converted to greyscale, passed through an edge detection filter, and finally converted to a true monochrome black/white bitmap. The end result of this preprocessing looks like this.

My goal is to trace a closed curve around the white pixels.

My naive approach was to make a list of all the white pixels in the image, then use a greedy travelling-salesman solver to draw a path around them. Which results in a path looking like this.

This works, there's a lot of artifacts but it works. But it's very slow. This relatively simple image takes a couple minutes to process. Any images more complex than this one need to be reduced in quality by a substantial amount before the time taken is reasonable. The resulting path is then almost unrecognizable.

When looking up algorithms for tracing images, I find a lot of things either built into image editing software (I need an algorithm that I can implement in Python), or overkill bitmap-to-vector algorithms designed to preserve color and line thickness.

All I need is a list of pixel values that trace out the lines.

There's got to be a better way than straight up solving travelling-salesman, right?

Depending on the accuracy of the trace you wish to achieve, there are a couple of ways you can do this:

The most straightforward way is to:

Convert the image to grayscale
Threshold, to achieve a completely binary image
Apply a high pass mask to get rid of the constant areas

For example $\begin{matrix} 1 & 1 & 1 \\ 1 & -9 & 1 \\ 1 & 1 & 1 \end{matrix}$

The result looks like this:

This technique is not without limitations, because as you can see, Pink Panther's whiskers have been drawn boldly and therefore have been "traced" and the highlight in his left eye has been interpreted as "background", producing an inner trace. To reduce these artifacts you can use a "wider" filter (for example, a $5 \times 5$ square of $1$s with a $-25$ in the middle. This produces a stronger wider trace line), morphological operations or adaptive application of the mask depending on features of the patch you are looking at.

The "path based" tracer you seem to suggest is how active contours work. In this case, you would initialise the active contour on the image rectangle and then iterate it towards the outer edges of the cartoon. Scikit-image has an implementation of active contours.

But based on the work you have done already, you could switch from a "Traveling salesman" route discovery to a convex hull of the points you collect.

The convex hull will trace an outline that encloses the set of points and give it a "cartoonistic" sketch view as well. Again, scipy has a convex hull implementation and so does Octave, from which I am posting this indicative output.

You will need to have the image package loaded for the following code to work and I am re-using the last image of the traced Pink Panther from above to extract the points, but that is just because I had it already, you can simply feed the convex hull the points you have extracted right after thresholding. It might be a lot of data indeed but convex hull is not solving a TSP and returns much quicker.

I = imread("thetracedimage.png");  % Load the image
Z = find(I);  % Find the outline (anything that is nonzero)
[Qy, Qx] = ind2sub(size(I),Z);  % Convert the index values returned by find, to image coordinates (in Qx, Qy) NOTICE ORDER OF X,Y!
H = convhull(Qx, -Qy);  % Obtain the convex hull
% H is the index of the existing points that make up the convex hull and we are assuming that these are connected by lines. To retrieve the actual points:
Fx=Qx(H);
Fy=Qy(H);
% Now plot everything. Notice the reversal of y purely because of the plotting coordinate system that puts 0,0 at the bottom left corner.
plot(Qx, -Qy,'.',Fx,-Fy);axis image;

And this results to:

Which of course is less accurate, but it depends on what you are trying to do and you seem to indicate that the approximate output of the TSP is sort of satisfactory (?).

Hope this helps.


Why motion cannot be estimated in DCT domain using block match? In all video coding standards, motion estimation is always done in pixel domain using block matching. I wonder why it cannot be estimated in DCT(Discrete cosine transform) domain using block matching directly? If two blocks are similar in value in DCT domain, doesn't it mean that they are similar to each other in pixel domain as well?

All the comments and suggestions are welcome!!!
Thank you!

Block matching will only work well if the pixel offset is (close to) a multiple of the block size, in both dimensional axises.  Odds against that might be (8*8-1):1 or (16*16-1):1.


Vein extraction with intersections I am not exactly working with blood vessels, but my problem is similar enough that I think someone who does work with blood vessels can answer this question for me.

I have several GPS points of a car driver going through a city. I would like to reconstruct the city map from these points.

The problem is that, unlike traditional vein extraction, there are intersections.

This is an example of the data I have: (reproduced as an image, but I have the actual coordinates.)

This is ONE example of what I would like it to become. (I have only segmented some of the roads.) I would be happy however if the algorithm classified each intersection as the start of a new road.

I can and I have used image processing to reproduce a contiguous map (using dilation and skeletonization). The process is not perfect, but if the algorithm works only with lines, I can produce the map in that form.

Note: (1) the map is much bigger than this, and it is unfeasible to do it by hand. (2) And no, I cannot use a real map.

Related posts:

Vein extraction from this image

Your problem is quite different than the vein extraction one. In fact, it is much mucg easier to solve with the "geometric" definition of the Hough Transform.

The Hough Transform is used exactly to detect lines in an image. It achieves this by integrating the brightness values of an image towards some angle $\theta$. So, for example, for angle $\theta = 0 ^ \circ$, the Hough Transform produces the sum of each row in the image, for angle $theta = 90 ^ \circ$, the Hough transform produces the sum of each column in the image...and so on.

The combined result of this is that points map to arcs, lines map to points and polygons map to specific configurations of points in the Hough Space.

The way to pick up a line in the image, in the Hough Space, is ultra simple. You just pick a maximum. But, that would give you an infinite line, not a line segment, which is what you have in your application.

The most straightforward way to pick up a line segment by the Hough Transform is to "mark" the pixels that are "responsible" for a particular accumulation. So, later on, when you pick up a maximum peak, you can "look" at the pixels that formed that sum and then connect them with a simple distance based rule. If distance is smaller than 4 pixels between any two pixels, then consider them as part of one line, otherwise create a new line.

For more information about all this, please see this link.

We now come to your case, you have two problems:

Detect line segments from GPS points
Connect line segments to form roads.

To solve #1, you can use the "geometric" definition of the Hough Transform WITHOUT having to work on the pixels of an image. Here is how this works:

Find the centroid of all GPS points. Call this the centre of your region.
Find the bounding box of all of your GPS points. These two steps define your "image space".
Run a Hough Transform with the center and over the bounding box of the points. As you are shooting "rays" down each different bearing, make sure that you also store the points that are responsible for the sum. The sum is simple: Shoot a ray down a direction, whenever it hits a point increase a counter and store its location.
Do a histogram of your Hough Space. It will help you in determining the threshold you have to apply so that you select the peaks that result in lines.
Apply the Threshold by selecting peaks that are above the thrshold value. You now have a bunch of infinite strong lines in your hands.
Look at the lists of points accumulated in each one of the directions and "integrate" them in lines if their distance is smaller than some threshold. You now have a bunch of line segments in your hands. These are your roads OR your road segments. The good thing about them is that they share a common point when they are connected (the junction). The bad thing about them is that they are unconnected when there is a "turn" or "curve".

You now have to solve problem number 2. Unfortunately, it is impossible to solve problem #2 without yet another application of a threshold. The philosophy here is simple:

Pick up a line, call it $a$.
Pick up all of the lines $a$ shares a common point with, call them $C$.
Calculate all angles between $a$ and the lines in the set $C$.
Merge two lines into one "road" if the angle is greater than some threshold. Obviously, if two lines form a $180^\circ$ angle, they are the same "road" but if they form a $90^\circ$ angle then they probably are not EXCEPT if set $C$ only has 1 element (1 line).

At the end of this process, you will have all connected segments that form roads.

It is not a difficult problem but it sure is a messy one.

Hope this helps.


How to Implement and Minimize Maximum Likelihood Expression in MATLAB I need to check if the estimation algorithm has converged or not. I am using the Maximum Likelihood estimation method. For convergence check, we see if the log-likelihood has reached its maximum value or not. But I am having difficulty in implementing the log-likelihood expression.

As an example, I am estimating the model parameters of a Moving Average model of order d =3 expressed in Eq(1). The known coefficients are h = [1    0.45   -0.2].The pdf is given in Eq(2) and the log-likelihood in Eq(3).

I have tried the following steps. For 9 different values of the channel coefficients, I calculate the likelihood expression. I then iterate 9 times to get 9 different values. After the vector of likelihood values are obtained, I should take the log of the values of likelihood and plot. However, I am stuck in the plot as I cannot understand how to show the maximum for the coefficients of the channel which is in a vector.

Question 1:  I am stuck in implementation in the step. How do I plot the curve and show at which value of $H$ the curve reaches its maximum value?

Question 2: I want to add noise of 5dB using the awgn() function as follows

 y=awgn(y,5,'measured');

But in that case I cannot understand what value I should take for the standard deviation of the noise, s?

 clear all

N=256; %Number of Samples to collect

H_est = [1.07648398119767,0.469220009437168,-0.0245459792881367;
1.01925235521823,0.449715366555421,-0.0287160909270392;
1.05903405523260,0.568022131010140,0.0261264158909992;
1.03693902020750,0.421708336037157,-0.102062737844706;
1.07969493608754,0.481505366768731,-0.0505480628598221;
1.05948950849694,0.450984209695615,-0.110305164694739;
0.933641686721610,0.310664477120195,-0.0803322291311877;
1.13980384345926,0.530800706049278,-0.0314494051827266;
0.949926376465499,0.498110774466918,0.0473466083388137];
 %Generating source information signal
z = rand(1,N);
h = [1    0.45   -0.2]; %true known channel impulse response

y = filter(h,1,z);
y=y+randn(1,N);
s=1; %Assume standard deviation s=1
tol = 1e-4;
maxiter = 9; %because there are 9 values of the channel coefficients
llh = -inf(1,maxiter);

 llh=zeros(1,9); %Place holder for likelihoods

for iter=1:9
%Calculate Likelihoods for each parameter value in the range

 H = H_est(iter,:);

     Hz=  filter(H,1,z);
        llh(iter) = exp(-sum((y-Hz).^2)/(2*s^2));  %Neglect the constant term (1/(sqrt(2*pi)*sigma))^N as it will pull %down the likelihood value to zero for

      %  increasing value of N
        if abs(llh(iter)-llh(iter-1)) < tol*abs(llh(iter-1)); break; end   % check likelihood for convergence
  end
        [maxL,index]=max(llh); %Select the parameter value with Maximum Likelihood
        display('Maximum Likelihood of H');
display(H_est(index,:));

plot(1:9,log(llh));

OK, let's see if I can answer this now.

The original log likelihood expression is
$$
L(\mathbf{y} | \mathbf{\Theta}) = -\frac{N}{2} \ln (2\pi \sigma^2_w) - \frac{1}{2\sigma^2_w} \sum_{n=0}^{N-1} ( y_n - \mathbf{h}^T \mathbf{z}_n)^2
$$
where $y_n, n=0,1,\ldots,N-1$ are the known noisy measurements, $\mathbf{h}$ are your (unknown) filter coefficients and $\mathbf{z}_n$ are the known random inputs to your filter.Here $\mathbf{\Theta}$ is just your unknown $\mathbf{h}$ coefficients if that is all that is to be estimated.

So, the way to calculate it is just (in R) as:

N <- 256
h_true <- c(1, 0.45, -0.2)
z <- rnorm(N,0,1)
hz <- filter(z,h_true)
sigma_w <- 1

y <- hz + rnorm(N,0,sigma_w)

llh <- -N/2*log(2*pi*sigma_w*sigma_w) - 1/(2*sigma_w*sigma_w)*sum((y-hz)^2)

Now, to your question:

  Question 1: I am stuck in implementation in the step. How do I plot the curve and show at which value of $H$ the curve reaches its maximum value?

Well, that's a little hard because your log likelihood is a multi-dimensional plot (four dimensional in your specific case -- one for each coefficient and one for the value of the log likelihood).

So you need to plot llh on three dimensions H_est(:,1),  H_est(:,2), and  H_est(:,3)  and then the value. You might want to just look at plot3(H_est(:,1), H_est(:,2), llh) (or something that approximate this).


Region growing algorithm I have been trying to come up with a region growing algorithm but I'm not sure that I fully understood the region growing segmentation method (for gray-scale images). Please correct me if I'm wrong:

I start from a seed point chosen by me (brightest value that fits the wanted region ,because the segmentation target is a girl's face).
From the location of that seed point I compare the adjacent points and if their value is between a given threshold and the seed point value then I add it to the seed point map(matrix).
I must try all the valid seed points for adjacent neighbors. The result is just plotting the seed point map (matrix).

Almost there:

Define a function growRegion(Image, seedX, seedY) that checks the pixel at location seedX, seedY and returns a list of the 4 or 8 pixels surrounding the seed that satisfy the criteria for inclusion. In your case, that is just one threshold so you could easily say, "if the pixel is adjacent to the seed pixel AND it is above the threshold then include it to the list of potential points".
Establish a list of points to be checked. For instance seedPoints=[(seedX,seedY)].
As long as seedPoints is NOT EMPTY

Pop a point (someSeedX, someSeedY) from seedPoints (remove it from the list)
Call the function returnPoints = growRegion(Image, someSeedX, someSeedY)
For each point in returnPoints mark a white pixel on some Mask image.
Extend seedPoints with the result of growRegion.

At this point, your Mask image contains all 4-8 connected pixels that satisfy the criteria for inlcudes. seedPoints fills up 4 (or 8) times the rate at which it empties so for large regions this might seem to run slowly.

Hope this helps.


Obtaining the size of objects from an image My objective is to be able to obtain the size, in pixels, of objects in images I'm receiving from a camera.

I've tried two approaches thus far, both working to some extent, and I'm wondering if there's any other approach you can suggest.

The first approach was using the CIE Delta E formulas. I would parse through the BGR 24 image data and then locate the appropriate pixels that form part of the object I wish to size. This approach works well when I know the object colors and positions beforehand, not so much when I don't. Another drawback of this approach is that calculating the Delta E for every pixel is extremely computationally expensive.

To circumvent my reliance on having to know the colors and positions beforehand, I used Connected-component Labeling to find the objects and size them.

Here's an example of what I'm doing using connected component labeling.

As you can see from the image above, I'm making four measurements, horizontal, vertical, left diagonal, and right diagonal. After I've sized the object I know where it is in the image, I then debayer the image into a BGR24 image and take several color samples. Debayering the image constitutes about 50% of the processing time, which is more than I prefer.

Now, my question is: What other approaches are there for finding the size in pixels of objects in an image?

I would use a split-and-merge algorithm for segmenting and then simply count the pixel per segment. Details are provided here

Two assumptions are neccessary:

1) you objects do not overlap.

2) the picture does not have any background or further irrelevant objects


Adding noise to an ECG signal I am doing a project on ECG arrythmia analysis using matlab.

I have designed notch filter for removing 50 Hz noise but don't know how to add a 50 Hz powerline interference noise to a clean ECG signal?
Also, I want to check whether noise is reduced in the filtered signal. Will Power spectral density using modified welch periodogram indicate whether noise is filtered or not?
How can I compare which wavelet (e.g. db6) is best suited for ECG analysis?

1) Create a 50 Hz sinusoid and then simply add it to your ECG signal.  You can control the power of the 50 Hz noise by multiplying the sinusoid by some gain factor (can be less than or more than 1) before you add it to the ECG.

2) I'm not familiar with the Welch periodogram, but if it displays the power spectral density then it should do fine.  I would just do an FFT myself.


Why do we care if a system is BIBO stable? I'm wondering what the actual application of this tidbit of knowledge is. I am currently a student in a Computer Engineering undergraduate program and I was wondering what the applicability of this is. Has there ever been a situation where a signal did come in bounded and leave unbounded, causing massive destruction?

The acronym BIBO, when used by university DSP professors, is merely a fancy (academic) term that means stable. So when when you encounter the term BIBO in the literature of DSP you can replace that term with the word "stable".

As for your final question, the answer is yes. Think of this: You have an audio signal that you'd like to apply to the input of a lowpass filter. And it's possible to design an IIR lowpass filter (a recursive filter, with feedback) using commercial software using filter coefficients that are represented by super-high precision floating point binary numbers. Then when you quantize your coefficients to 8-bit values in order to implement your lowpass filter using an 8-bit microcontroller, it's possible for the "imprecise" 8-bit filter coefficients to cause the IIR filter to become unstable.

That scenario, described in every DSP textbook, means each input sample applied to your filter will cause the filter's successive output samples to uncontrollably grow in value. And when fixed-bit-width binary numbers grow in value they will eventually experience overflow error.  The bottom line here is: When an information-carrying discrete signal is applied to an unstable digital filter the output quickly becomes random noise. To quote my young granddaughter, "This is a bad thing."


camera calibration: why chessboards? Chessboards are very frequently used for camera calibration, this website also talks a bit about why: wiki. One of the reasons is, that chessboards define many natural interest points.
Now my (very simple?) question: Why don't we just use a white board with black dots on it instead of the checkerboard pattern where we first need to detect the corner points? I believe the black dots would be even easier to find...

This is mostly a matter of history.

Checkerboard patterns have proven in the past to:

yield calibration accurate results
be easy and robust to implement.

This was confirmed (for example) by a benchmark performed by NASA long time ago, where they were assessing visual solutions to control robotic arms.

Why are checkerboards robust and accurate? Because they provide features that can ve computed as intersection of lines.
Linear features were easy to detect accurately (using Canny-Deriche, SUSAN, Hough transform...), while most point detectors were computed with pixel accuracy.
Furthermore, checkerboards do provide line intersections, which are:

robust: if a part of the calibration pattern is occluded, you can still get the coordinates of the intersection points thanks to the line equation
accurate: line equations can be interpolated from sample points, and this is the case for their intersection too. Thus, you naturally get subpixel precision. This holds also for the old-fashioned pipelines where the corners were actually clicked by somebody

By robust to occlusion, I mean that you just have to continue the parts of the lines that are visible to get the intersection coordinates.
With dots this is more hazardous, because you have less clues to infer the overall line direction. Also, when some dots are occluded, it's very difficult to know which ones are indeed occluded. Instead, using checkerboards of known size (e.g. $8 \times 8$ squares), then you can infer which tile exactly is hidden as long as you see the borders of the checkerboard.

A possible advantage of checkerboard patterns is that since you have line featurs, then it may be easier to infer distortion coefficient values, which is harder to do if you need to "connect the dots" of aligned points.

Modern point detectors (starting from Harris corners) perform interpolation to reach subpixel accuracy in point detection, and since checkerboards provide a lot of corners they continued to be used, just switching from line primitives to corner primitives.


Why can we know the best precoding direction after we do the SVD to the channel? Why can we know the best precoding direction after we do the SVD to the channel?for example

$\mathbf y=\mathbf H\mathbf s+\mathbf n$,now we do the SVD to the $\mathbf H$

If the $\sigma _1$ is the biggest singular value,after we normalize $V_1$,why can we say the $V_1$ is the best precoding direction?Is there any theorem can prove this?

First you have to know that the vector in $U$ and the vector in $V$ are orthonomal,that is ,$\vec u \vec u^H=I$,and $\vec v \vec v^H=I$.

So now $H=U\Sigma  V^H$,and $y=Hs+n=(U\Sigma  V^H)s+n$

and $U^H y=\Sigma  V^Hs+n$,see the $\Sigma  V^Hs$,isn't it like " channel gain $\times$ beamforming $\times$ signal $s$ "?

Now you can also find that according to the rule of multiplication of matrix,the first element of $\Sigma$,which is the $\sigma_1$,will multiple every element of $v_1^H$,that is why we say the $v^H_1$ is the best precoding ,because its elements will multiply the biggest singular value ,$\sigma_1$, every time.


How to Determine Covariance Matrix $Q$ and $R$ in Kalman Filter I am implementing getting orientation from smartphone. I want to use Kalman filter and should determine process noise covariance matrix $Q$ and measurement noise covariance matrix $R$. (newbie to Kalman filter)

I don't have any idea how to determine $Q$. What I think about $R$ is as follows:

state vector : quaternion from (accelerometer + gyroscope)

(1) My phone is stand still. I get covaraince matrix from Matlab

1.0e-04 *

0.0000    0.0005    0.0035   -0.0000
0.0005    0.0063    0.0411   -0.0002
0.0035    0.0411    0.2881   -0.0014

-0.0000   -0.0002   -0.0014    0.0000

(2) My phone had been moved for 5 seconds.

covariance matrix is

0.0417   -0.0533   -0.0008   -0.0014
-0.0533    0.0784    0.0015    0.0018
-0.0008    0.0015    0.0001    0.0001
-0.0014    0.0018    0.0001    0.0001

Is there anyone to help?

(added)

Details are omitted.

case 1: Kalman Filter

The row data from my phone is $p, q, r$ (angular velocity). I omit the conversion equation between angular velocity and quaternion.
\begin{align*}
x_{k+1} &= Ax_k+w_k \\
z_k &= Hx_k + \nu_k \\
Q &: \text{ covariance matrix for }w_k\\
R &: \text{ covariance matrix for }\nu_k
\end{align*}

$$ \begin{bmatrix}\dot q_1\\\dot q_2\\\dot q_3\\\dot q_4\end{bmatrix}=\frac 12 \begin{bmatrix} 0&-p&-q&-r\\p&0&r&-q\\q&-r&0&p\\r&q&-p&0 \end{bmatrix}
\begin{bmatrix} q_1\\ q_2\\ q_3\\ q_4\end{bmatrix}$$

$$ \underbrace{\begin{bmatrix} q_1\\ q_2\\ q_3\\ q_4\end{bmatrix}_{k+1}}_{x_{k+1}}=\underbrace{\left( I + \Delta t\cdot \frac 12 \begin{bmatrix} 0&-p&-q&-r\\p&0&r&-q\\q&-r&0&p\\r&q&-p&0 \end{bmatrix}\right)}_{A}
\underbrace{\begin{bmatrix} q_1\\ q_2\\ q_3\\ q_4\end{bmatrix}_k}_{x_k} $$
$$ H=I$$
My guess for covariance matrix is as follows: (but I don' know how to infer..)
$$ Q = 0.001I, \quad R=10I.$$

case 2: Extended Kalman Filter

\begin{align*}x_{k+1} &= f(x_k) + w_k \\
z_k &= h(x_k) + \nu_k\\
Q &: \text{ covariance matrix for }w_k\\
R &: \text{ covariance matrix for }\nu_k
\end{align*}
$$ A = \left.\frac{\partial f}{\partial x}\right|_{x_k} ,\quad H = \left.\frac{\partial h}{\partial x}\right|_{x_k} $$
\begin{align*} \begin{bmatrix} \dot \phi\\ \dot \theta\\ \dot \varphi \end{bmatrix}&=
\begin{bmatrix} 1&\sin\phi\tan\theta & \cos\phi\tan\theta \\
0&\cos\phi & -\sin\phi \\
0&\sin\phi\sec\theta & \cos\phi\sec\theta \end{bmatrix} \begin{bmatrix}p\\ q\\  r \end{bmatrix} \\
&= \begin{bmatrix} p+q\sin\phi\tan\theta+r\cos\phi\tan\theta \\
q\cos\phi-r\sin\phi \\
q\sin\phi\sec\theta + r\cos\phi\sec\theta \end{bmatrix} \\
&= f(x) + w
\end{align*}
$$z = \begin{bmatrix}1&0&0\\0&1&0\end{bmatrix}\begin{bmatrix}\phi\\\theta\\\varphi\end{bmatrix}+\nu = Hx + \nu $$
$$A = \begin{bmatrix} \frac{\partial f_1}{\partial\phi} & \frac{\partial f_1}{\partial\theta} & \frac{\partial f_1}{\partial\varphi} \\
\frac{\partial f_2}{\partial\phi} & \frac{\partial f_2}{\partial\theta} & \frac{\partial f_2}{\partial\varphi} \\
\frac{\partial f_3}{\partial\phi} & \frac{\partial f_3}{\partial\theta} & \frac{\partial f_3}{\partial\varphi}
\end{bmatrix} $$

(I emphasize that details are omitted.) In this case, also I don't know how to infere $Q, R$.

R depends on the sensor sensitivity. If this is a real world problem this can be obtained from the manufacturer. If not use the identity matrix multiplied by a scalar that is less than 1.
Q is the covariance of the process noise. Again if this is a real world problem this can be obtained in the noise level in the states of the system at steady state. if not you can assume Q is zero matrix. A non zero Q helps achieve good convergency characteristics as explained here.


How to find poles of transfer function by looking at the step response? How to find poles of transfer function by looking at the step response?

Given a step response graph like such:

How would I find the sketch for its poles on the complex plane? The only thing I can gather are the rise time, settling time, overshoot, etc. But how do I derive the poles from that?

There are two things to determine: the frequency $f_n$ (or angular frequency $\omega_n$) and the $Q$.  Approximately (for a high $Q$ system like the one pictured):

1) period $T$ is the time between successive peaks. $f_n=1/T$, and $\omega_n=2 \pi f_n$.

2) ratio of successive positive (or negative) peaks (measured relative to steady-state) $ 1- \delta = e^{-2 \pi /(2Q)} $, or, for high $Q$, $ Q \approx \frac{\pi}{\delta}$.  [Example: If $\delta=0.1$ (10% amplitude decay per cycle), $Q \approx 30$.]

The poles are at $ \omega_n \left[-\frac{1}{2Q} \pm i \sqrt{1-\frac{1}{(2Q)^2}} \right] $, or in radial coordinates: $(\omega_n, \pi \pm \tan^{-1}{\sqrt{4Q^2-1}})$.

[There's also a pole at zero for the step response. Per bobuhito's comment, this pole is not part of the transfer function.]


Find the transfer function of the difference between IIR and FIR filter I am using a filter following equations from papers.
It is basically the difference between exponential moving average and simple moving average.

$$L[n]=\frac{1}{\alpha_L}f[n]+\left(1-\frac{1}{\alpha_L}\right)L[n-1]$$

$$ S[n]=\frac{1}{\alpha_S}\sum_{k=0}^{\alpha_S}f[n-k]$$

$$ g[n]=S[n]-L[n]$$

$f[n]$ - original signal
$g[n]$ - filtered signal
$\alpha$ - window size

I can find their transfer function separately for $L$ and $S$, but I have no idea how to get transfer function of $g[n]$. Anyone please guide me, I would really appreciate it.

Because you are using a linear combination of the two filters, you can find the resulting transfer function by directly subtracting them from each other. If you want to do this, then you have to make sure that they have the same denominator, such that the total transfer function is one quotient of polynomials. If you do not know how to do this, you could check this example.


Why Do I Get This Crackling Noise on Zeroing out the High Frequencies? I recently started playing with the Fourier transform (after spending a few weeks learning about the mathematics behind it). I decided to try hacking together a low-pass filter on the following sound bite:

In particular, I took the Fourier transform, zeroed out the highest 1/2 of the frequencies, and then took the inverse Fourier transform. This is what I got

Why is there that crackling noise?

Two potential issues:

Doing filtering in the frequency domain (using an FFT) requires overlap-add, overlap-save, or a related algorithm. This is caused by the difference between linear and circular convolution. Otherwise you get time-domain aliasing
Some of the noise sounds like simple clipping. Filtering can actually increase the time domain amplitude for some samples and if that exceeds the available range it either clips or wraps around.


Aliasing in digital communication system In a digital communication system, is it possible to counteract the aliasing problem by properly tuning the constellation order $M$ and/or the carrier frequency?

Here is my train of thought, please correct me if I'm wrong:

The sampling rate needs to be 2 times the signal bandwidth. This is because we need to cover all the frequencies in which the signal operates. If the signal breaks the Nyquist criterion we'll only get a part of the signal. So if I look at a fix sampling rate and try to tune the constellation order $M$ up and down it won't help, because the signal will still have the same bandwidth. The only thing that will change is the transfer rate $r_b = \log_2(M)/T_s$, where $T_s$ is the sampling rate. Also tinkering around with carrier frequency will only move the signal up or down the frequency. So neither of these solutions will counteract the aliasing problem. The only solution is to somehow reduce the signal bandwidth. I'm thinking of using an RRC filter on the transmitter (to shrink the bandwidth) and then using another RRC filter on the receiver to counteract the ISI.

If the signal breaks the Nyquist criterion we'll only get a part of the signal.

This needs to be clarified: the frequencies above Nyquist will fold back into the Nyquist band, so it's not just that you lose part of the signal; it's that the signal is corrupted (usually) beyond hope of reconstruction.

As an example, consider a 4 kHz real signal sampled at 4 kHz -- the frequencies between 2 and 4 kHz will fold back into the 0 to 2 kHz range; not a single "part" of the signal will survive.

You're correct that this can't be solved by modifying the constellation order $M$ or the carrier frequency. As you correctly assume, you have two options. You can reduce the bandwidth through pulse shaping; if the symbol rate is $R$, then the minimum bandwidth $B=R/2$ is achieved with sinc pulses. The second option is to reduce the rate.

(As a side note, recently there have been proposals to increase the rate without increasing the bandwidth (search for "Faster than Nyquist signaling"), but I think these ideas are outside the scope of your studies at the moment.)


Formant Analysis Frequency - What makes a peak value? from spectrograms, how can you tell which peak is really a formant. Like why is the pink arrow the first formant. And in the spectrogram below, which has fewer samples, what makes a value the peak value?

The spectra on the left are two line spectra of two different periodic functions, the top with fundamental frequency of 100 Hz and the bottom has $f_0$ of 200 Hz.  While they are not perfectly flat spectral lines, they are broadbanded with no apparent resonance.

It's a gross approximation, but I believe these two broadbanded periodic signals are meant to represent what comes out of one's voicebox in the throat.  It's the same broadbanded, except the $f_0$ depends on the pitch that one intends to sing (or speak) at.  If it were A above middle C (normally called "A4" or MIDI note 69) then $f_0$=440 Hz.

Now that broadbanded (same spectral envelope, but different $f_0$) driving signal on the left gets passed through an acoustic filter in the middle which is not dependent on $f_0$ but is dependent on the shape of the mouth, position of tongue, amount of nasal coupling, all which is dependent on the vowel one is singing.  In the frequency domain, the spectrum of the driving signal is multiplied by the frequency response to become the spectrum of the output.

The cause of the formant frequencies is due to the acoustic filtering of the vocal cavity that is one's head, mouth, lips, and nasel coupling.

A few of us have developed a pitch shifting algorithm that shifts pitch, but does not shift the formant frequencies.  I did this a quarter century ago and wrote a short paper about it.  You can sorta think of your voicebox as emitting a quasiperiodic string of impulses (or approximations to impulses) from the glottal closure and opening (with fundamental frequency $f_0$) and those impulses being convolved with the impulse response of the formant filter.  If the pitch shifting does not stretch or scrunch that impulse response, the formants remain unchanged.

But simple speeding up or slowing down the recording does scrunch or stretch the impulse response, which causes the formants to move.  That causes someone to sound like Alvin the Chipmunk (pitching up) or Darth Vader or Satan from hell (pitching down).


How to reverse Mel scale The conventional MEL scale converts the frequency scale using the following formula:

$$M_f = 1127.01048 \cdot \log\left(\frac{f}{700} +1\right)$$

I want to use a reverse scale where resolution is better at higher frequencies only up to $10\textrm{ kHz}$.

In my opinion it will be a s-shaped curve.
Any suggestions?

The formulation for MEL frequency scale is logarithm scale based, it clubs(put together) the higher frequencies more closely(so less resolution), because logarithmic functions vary slowly for larger numbers(see the logarihmic functions plot).

So some function of exponential form will definitely work out if one wants to give more resolution to higher frequencies.

so based on particular requirements one can define exponential scale for effect that looks like reverse of MEL scale.

Following is the filter bank I created using exponential function, similar to the formulation of MEL scale.

$fnew$ $=$ $exp(f / 7000)$.

I think this is what you are looking for, more resolution for higher frequencies as compared to lower frequencies.


How come FFT window give half size frequency bands when filter can seems to have an infinite precision? I don't really understand maths behind filtering but use it lot in music production and i'm a total newbie to DSP.

Usually parametric filters let you filter frequencies with a precision that seems almost infinite. But when using a FFT or DFT the taken window size gives n / 2 + 1 frequency bands.

Are the two related ? if so how can a filter get so much precision when Fourier Transform doesn't?

That may seem a dumb question but when I look at spectrograph in sound softwares they always looks much more precise than the buffer I've choosed (usually 256 samples).

The DFT is a tool to represent the spectrum of any signal. Usually for an N point DFT only N/2 output coefficients are required because for a real signal the DFT is symmetric about the centre.

No filter has infinite precision, neither infinite bandwidth nor numeric representation. A filter with enough precision can be chosen for an application to seem smooth enough.

The DFT though doesn't have anything to do with a filter. It simply represents the spectrum of the signal at it's input.


Speech comparison algorithm for rating on similarities I am trying to compare 2 speech samples and rate them on similarities. Think of someone trying to repeat a phrase, and then comparing those 2 audio files.

I started by implementing the MFCC (http://en.wikipedia.org/wiki/Mel-frequency_cepstrum) algorithm. I calculate the MFCCs of both audio samples, which gives me roughly 500 frames of audio (at 10ms each, with like 30% overlapping of the previous) having 14 or so MFCC coefficients. So a 500x14 matrix for each audio signal.

Then I do the naive approach of simply differencing the matrices. This does not give very promising results. Half of the time when I compare completely different audio samples (where different phrases are spoken), I get less difference than comparing the audio where I try to repeat the same phrase! This is clearly backwards and can't give me a good scoring algorithm.

How can I improve on this? I thought MFCCs were a really important part of speech processing, though clearly I need to do more with it.

First, you will have to correct for differences in timing. For example, if one utterance is "--heeelloooo---" and the other "hellooooooo----" (- representing silence), a direct pairwise comparison of MFCC frames will show differences simply because the two samples are not aligned. You can use Dynamic Time Warping to find the best alignment between the two sequence of feature vectors - and compute the corresponding distance.

The second problem is that if the two recordings are not from the same speaker, you will have to compensate for the differences in timbre. The MFCC of a female saying "aaa" are not the same as the MFCC of a male speaker saying the same phoneme! A relatively simple model to account for variations in voice timbre is to assume that there exists a linear transform $\Gamma$ which "maps" the MFCCs of one speaker onto the MFCCs of another speaker (to be fair, only one small subset of these transform accurately model how changing parameters such as age, gender, etc "shifts" the MFCC). Once two recordings have been aligned, even roughly, you can use a least square procedure to estimate $\Gamma$. This procedure is known as speaker normalization or speaker adaptation.

Your comparison procedure will thus consists of the following steps. $A$ and $B$ are your original MFCC sequences.

Align the two utterances using DTW. This yields $A'$, a matrix with the observations from $A$ warped/shifted to be aligned with the observations in $B$. You can stop here if $A'$ and $B$ are known to be from the same speaker.
Estimate the transform $\Gamma$ which minimizes the difference between $\Gamma A'$ and $B$.
Use the distance between $\Gamma A'$ and $B$ as your metric.

A last thing that comes to my mind is that you have to discard the first MFCC coefficient (which roughly expresses the signal loudness) to improve the capacity of your system to match utterances pronounced at a different volume/recording level.


How can I use a trained model to generate related images? Suppose I have a data set of a large number of related images and I extract HOG features and train it using SVM Classifier or KNN Classifier. How can I use this model to generate a new image which is related to the images in the dataset (rather than testing on a new image which we do traditionally)?

There are multiple aspects which render your approach infeasible:

SVM is a learning technique where high-dimensional data is mapped to lower dimension (a couple of hyperplanes, decision boundaries, which pass through support vectors). It doesn't have an inverse mapping: It is not super trivial to go from a point within a decision space back to the space of actual data. In other words, it's not generative. If your learning method supports going back from low-dimensional embedding to the higher dimensions space (such as GPLVM), your task could be easier.
HOG is not exactly reversible. Please check Hoggles, where the authors reveal what HOG features tend to understand, via formulating the visualization as in inverse problem. Therefore, it might not be able to synthesize consistent colors, or even detailed features. The synthesized image might end up looking just like a Hoggle image.

I would suggest the Generative Adverserial Networks for this task, which have been quite successful in image synthesis:

Here are some papers about it:

https://arxiv.org/pdf/1612.03242v1.pdf
https://arxiv.org/pdf/1406.2661v1.pdf

Here are some open source software:

https://github.com/junyanz/iGAN
https://github.com/carpedm20/DCGAN-tensorflow
https://github.com/jacobgil/keras-dcgan


System invertability Can someone please explain how to determine if a system is invertible or not. The particular systems that I am looking at are:

$y(t) = \cos \left( x(t) \right)$ - I believe this is non-invertable because of the periodic function

$y[n] = nx[n]$ - I think this is invertable and the inverse is $x[n]/n$

$y[n] = x[4n + 1]$ - Is the inverse $x[1/4n - 1]$?

$y[n] = x[n - 2] - 2x[n - 8]$ - I do not know even how to begin

I am really confused with systems... did not understand a single thing in the lecture today, and now reading the book ("Signals and Systems", Oppenheim) I am getting frustrated. Is there an algebraic procedure for determining the inverse of a system?

Please help.

$$ y[n] = x[n-2] - 2 x[x-8] $$

that's actually a linear and time-invariant system.  you invert it by computing the Z transform

$$ Y(z) = z^{-2}X(z) - 2 z^{-8}X(z) $$

solve for $X(z)$, then inverse Z transform.  you might find out that, to invert this system, you will have to look at future samples, so the inverse would not be causal.


How Could I extract the signal in the below Image This is my first post so sorry if I make any mistake.

The entire signal length in the screen shot below is 10 second.
I ran a first order low pass filter (digital) at 10Hz cut-off frequency, and I obtained the orange/yellow line in the picture. Please ignore the red and the blue line for now.

I wonder if there is an algorithm that I can extract all of the lowest values of the green plot in real-time and reconstruct/interpolate data along the way like the black line that is shown below?

The green signal (50Hz) is the DFT result of a raw 1kHz sine data that is sampled at 10kHz. I really hope that we can find a light-weight algorithm for the filter, since the entire code has to run on a microcontroller.

I would love to hear your suggestions. And let me know if anything is needed for clarification.

Normally you could do some valley detection. Scipy does this (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html) but this requires knowing a chunk of the signal in advance so you couldn't do it real time.

For real time maybe you could look at the derivative of the signal and notice when it goes from negative to positive, it would find valleys although it may not be as smooth as your black line. Maybe also check out this: https://stackoverflow.com/questions/22583391/peak-signal-detection-in-realtime-timeseries-data


2D Digital Signal Processing - Applying High Pass on and Image I have been doing my research on 2D digital signal processing. What I am doing is about noise removal of digital image:

Type of noise that I want to remove from an images is like the Gaussian noise and 'salt and pepper'. For your information, I intend to use the butterworth lowpass filter to remove the noise. This is my code that I used in matlab.

Function Highpass

function [ out ] = butterhp( im,d,n )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
h=size(im,1);
w=size(im,2);
[x,y]=meshgrid(-floor(w/2):floor(w-1)/2); -floor(h/2):floor(h-1)/2;
out= 1./(1.+(d./(x.^2+y.^2).^0.5).^(2*n));
 end

Function Lowpass

function [ out] = butterlp( im,d,n)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
out=1-butterhp(im,d,n);
end

matlab code

a = imread('penang.jpg')
f = rgb2gray(a);
k = imnoise(f,'gaussian')
lp=butterlp(k,15,1);
af=fftshift(fft2(k));
aflp=af.*lp;
aflpi=ifft2(aflp)
noise = uint8(real(aflpi));
imshow(noise)

This is the image that i used.

And running the code, this is what I get.

So,my problem is when I try to filter it to make it looks more clear by adjusting the number of 'n-th order' and 'cutoff frequency', it seems it does not work. It still gives me blur picture just like that.

First things first: If you use fftshift, you have to undo it before doing the ifft. Try

af=fftshift(fft2(k));
aflp=af.*lp;
aflpi=ifft2(ifftshift(aflp))

The result should improve, especially this checkerboard-overlay should go away.


Soft bits in large constellations I am working on a soft iterative receiver for an M-ary modulation.
Since my constellations, M=2⁷ — M=2¹², are rather big, the log likelihood ratios also become large.
Is there a way to compress these values to a smaller range, e.g. from [-100, 100] to [-10,10]?
Momentarily, I am rescaling the values to fit them in my desired range. However, the actual confidence levels don't make sense then anymore.
== Edit ==
To calculate the soft-outputs of my decoder, I first convert the soft-inputs into probabilities by,
$$
P = \frac{e^{LLR_{in}}}{1 + e^{LLR_{in}}}
$$
As can be seen from this equation, large LLRs will be approximately 1. However if this happens, my decoder takes wrong decisions and calculates wrong soft-outputs. Therefore I am searching for a way to deacrease the magnitudes of the input LLRs.
The MATLAB rescale() function isn't the optimal solution, since it does introduce errors with the iterations.

It depends on your decoder; for example, in message-passing decoders, most commonly, a prefactor doesn't affect the results at all, so you can scale to your heart's delight.
It's also common to just clip – whether a bit has a LLR of -100 or -90 doesn't really matter, does it.

However, the actual confidence levels don't make sense then anymore.

Well, you'd have to transform them, too. Notice that scaling LLRs makes sense in a Gaussian noise setting, as LLR involves taking the logarithm of a $ e^{-a|x-x_i|^2}$ function, so you're left with something rather sensible (scaling the LLR just scales $a$, which is common to all samples and constellation points). I don't know how your noise looks or how you arrive at confidence levels, but try to apply the logarithm to a ratio of confidences and see whether scaling the LLR preserves sensibility.


Is real-time signal processing possible under Windows? This question is related to on How to decide what platform is best to implement real-time audio processing on?, but focusing on Windows.

I'd like to be able to take input from a data acquisition board and process it (either through Matlab or a custom program), and use the result of that program to drive an external device that is attached through USB.  I know that getting hard real-time responses with Windows is nearly impossible (as there is a "wiggle" (technical term) of 3-15ms).

What is the most effective way to decrease this wiggle time?  Is this possible?  Should I be looking for an embedded solution straight away?

if you need to keep jitter (assuming that's what you mean by "wiggle" :) down to <1ms I'd say forget Windows.  Does the absolute latency from input to output matter, or is it just the jitter?

You might manage it some of the time, but guaranteeing it is unlikely.  What's the penalty for failing to meet the deadline?  Will it blow up speakers?  Drive a robot arm too hard against an end-stop?  Or something more benign?

The fact that your output is USB is also not likely to help in the jitter stakes as there's lots of extra OS interaction there, and the topology of the bus will also hinder your potential latency - if your output is at the end of a string of a few hubs, for example.


Understanding the frequency domain I'm trying to understand what frequency domain is. I found general explanations on the Internet, for example:

frequency-domain graph shows how much of the signal lies within each given frequency band over a range of frequencies
Fourier transform refers to both the frequency domain representation and the mathematical operation that associates the frequency domain representation to a function of time.
Fourier series transform a signal from time domain to frequency domain.

But I could not find an example which shows how I can obtain the frequency domain graph using the formula of a signal.
The graphs I found are sometimes discrete like this one:

and sometimes continuous like this one:

that seems confusing.
Suppose I have a signal $x:\Bbb R\to \Bbb R$ defined by the formula $x(t) = \cos(6\pi t)e^{-\pi t^2}$.
It's clear how to plot its time domain graph, But how can I find its frequency domain function or graph?
Should I apply fourier transform on it?:
$$\widehat{x}(\tau)=\int_{-\infty}^\infty \cos(6\pi t)e^{-\pi t^2}e^{-2\pi t \tau i}dt$$
Does plotting this new function give the frequency domain graph? (but it is complex valued, how can it be plotted?)
Or should I find the Fourier series of $x(t)$ and plot the series coefficients discretely?
So my only question is: How can I mathematically obtain the the frequency domain function (and then plot it to get the frequency domain graph) using the formula of $x(t)$?

In simple terms, Image shown here speaks for itself.

Before speaking about Fourier Transform black magic lets understand idea behind it. Work of the Mathematician Joseph Fourier demonstrated any arbitrary periodic (this is important) signal can be decomposed into bunch of sine-waves at with their corresponding amplitude and relative phases.

Essentially, he meant basic component to build any signal is a sine-wave. In other words you can build any arbitrary signal exactly if you have information about how many different frequencies are present, amplitude of the each sine wave, and there relative phases. This job can be done by inverse Fourier transform. And as should be clear by now, Fourier Transform does the reverse give it composed time domain signal and it will tell you sine-waves amplitude for corresponding frequencies and there relative phases. Hence, the output is a complex. Which makes sense in polar co-ordinates, which gives you amplitude and phase.

E.G. When you strum a guitar strings, it vibrates at multiple frequencies with different amplitude which makes up a beautiful sound, otherwise you know how a single frequency sounds(watch this :https://www.youtube.com/watch?v=qNf9nzvnd1k). one with the highest amplitude is known as tone.

So, as you can see in the figure attached here, Fourier transform tells you exactly the same information.

  Suppose I have a signal $x:\Bbb R\to \Bbb R$ defined by the formula
  $x(t) = \cos(6\pi x)e^{-\pi x^2}$.

For to obtain Fourier transform of the function you specified it should be function of time (t). So, I will just re-write your function above slightly differently replace x with t, since x(t).

$x(t) = \cos(6\pi t)e^{-\pi t^2}$.

Also, Take a look at this,
http://www.thefouriertransform.com/pairs/rightSidedSinusoids.php

Hope this helps.


Inconsistency Between Analysis and Simulation Results of Fourier Transform of a Sinc Function First of all, hello.
This question is about a problem that I’ve faced during an attempt to obtain both time and frequency responses of a sinc function in MATLAB. The problem is an inconsistency between my calculations that I’ve done on the paper and simulation results that I’ve observed in MATLAB.
I’ve picked up $sinc(10 \cdot t)$ function and at first, I’ve taken its Fourier transform on paper. Afterwards, in order to have an experience in coding, I implemented this operation into MATLAB and plotted the functions both in the time and frequency domains.
However, in the amplitude spectrum in the frequency domain, the amplitude of the resultant rectangular pulse didn’t match with my results on the paper. According to my analysis results, the amplitude of the rectangular pulse is 0.1 but in MATLAB, I’ve got 10 instead of it.
My MATLAB code for simulation is given below:
t = -10: 0.01: 10;
Ts = 0.01;  % Sampling period in seconds
fs = 100;  % Sampling frequency in Hertz
function_1 = sinc(10 * t);
number_of_samples_1 = length(function_1);
number_of_samples_2 = power(2, nextpow2(number_of_samples_1));
f = linspace(-fs / 2, fs / 2, number_of_samples_2);
function_2 = fft(function_1, number_of_samples_2);
function_3 = fftshift(function_2);
function_4 = abs(function_3);
subplot(2, 1, 1);
plot(t, function_1, "r");
xlabel("Time (s)");
ylabel("Function-1 (t)");
title("Sinc(10.t) Function in the Time Domain");
subplot(2, 1, 2);
plot(f, function_4, "r");
xlabel("Frequency (Hz)");
ylabel("|Function-4 (f)|");
title("Amplitude Spectrum of Sinc(10.t) Function in the Frequency Domain");

My question is this: Why has this result occurred in contrast to my expectations and analysis?
I will be happy if you give feedback to my problem and question.
Thank you.

I have recently answered a related question.
What happens is that MATLAB fft is calculating a sum in the sampled points, the integral would be approximated with a (Riemann) sum.
$$ \int_{0}^{T} f(x)dx \approx \sum_{k=0}^{N-1} x(k \Delta t) \Delta t$$
The $\Delta t=0.01$ factor matches the two results.


Does a matched filter give a gain in SNR? I am trying to know if filtering a sequence with a matched filters gives a gain equal to the sequence length, he are the calculus, is it correct ?

In my model $s$ is a deterministic complex sequence (lets say RRC pulse shape samples for a QAM modulation), $v$ is a complex white circular gaussian noise.

If it is correct does it mean that there is an "inband" filter which can enhance $SNR$ ?!

\begin{align}
s &= [ s_0,\dots,s_{ N-1} ]^T\\
v &= [ v_0,\dots,v_{ N-1} ]^T \\
x &= s + v \\
R_v &= \mathbb{ E } \left[ v v ^ H \right] = \sigma ^ 2 I \\
\textrm{SNR}_1 &= \frac{ \mathbb{ E } \left[ \lvert s \rvert ^ 2   \right] } { \mathbb{ E } \left[ \lvert v \rvert ^ 2 \right] } = \frac{ \mathbb{ E } \left[ s ^ H s  \right] } { \mathbb{ E } \left[ v ^ H v \right]} = \frac{ s ^ H s } { \mathbb{ E } \left[ v ^ H v \right] }  = \frac{ s ^ H s } { \mathbb{ E } \left[ \Sigma \lvert v_i \rvert ^ 2 \right] } = \frac{ s ^ H s }{ N \sigma ^ 2 }\\
h &= s \\
y &= h ^ H x = h ^ H s + h ^ H v = s ^ H s + s ^ H v  \\
\textrm{SNR}_2 &= \frac{ \mathbb{ E } \left[ \lvert h ^ H s \rvert ^ 2   \right] } { \mathbb{ E } \left[ \lvert h ^ H  v \rvert ^ 2 \right] }
  = \frac{ h ^ H s s ^ H h }{ \mathbb{ E } \left[ h ^ H v v ^ H h  \right] }
     = \frac{ h ^ H s s ^ H h }{ h ^ H R_v h } = \frac{ s ^ H s s ^ H s }{ s ^ H \sigma ^ 2 I s } = \frac{ s ^ H s }{ \sigma ^ 2 } \\
     \frac{\textrm{SNR}_2 }{ \textrm{SNR}_1 } &= N
\end{align}

Output SNR for matched filter can be analysed in time domain, which is as explained below,

Consider a filter with Impulse response  $ \enspace h(t)$  and input $\,s(t)+w(t)$, where $\,s(t)$ represents the required signal and $\,w(t)$ represents AWGN noise.

Output signal is given by  $$y(t) = s(t) * h(t) + w(t) * h(t) $$

Writing signal component of output in inverse fourier form   $$ y'(t) = \int_{-\infty}^{\infty}  S(f) H(f)  \ e^{j2pit}\ df $$ where $S(f)$ and $H(f)$ are spectrums of input signal $s(t)$ and Frequency response of filter.

Output SNR of filter which is sampled at time instant $\,t=T $ is output power with signal component by noise power which can be return as

$$SNR  =  \frac{\displaystyle\left |\int_{-\infty}^{\infty}  S(f) H(f)  \ e^{j2piT}\ df \right|^2}{ \displaystyle\left|\int_{-\infty}^{\infty}  \frac{N_o}{2} |H(f)|^2  df \right|^2}$$

Note that numerator term is evaluated at time t = T, which is the sampling instant for Matched filter.

To maximise this SNR , we need to cancel the phase term .
if we substitute $$H(f) = S^*(f)  e^{-j2pif} $$
which is nothing but $$h(t) = s^*(T-t) $$

and simplify the above SNR expression, we will get
$$ SNR = \frac{\displaystyle\int_{-\infty}^{\infty} \left|S(f)\right|^2 df}{N_o/2} $$
which is nothing but ration of energy of signal $s(t)$ to value of power spectral density of AWGN noise.

In this time domain analysis you can see that there is no filter effect , and it also tells us what should be the filter impulse response to maximise the output SNR.


How is angle of arrival estimated? I have simulated a phased array antenna. I have also simulated a point source. I don't get how to calculate pseudo spectrum vector without using the position of the source. Algorithms such as Capon and Bartlett use the position of the source for calculate the auto- correlation matrix which is further used to calculate pseudo spectrum. How is the angle of arrival being estimated if the source position is required to run the algorithm in the first place?

Here is the code:

M=6;

sig2=.1;

th1=-3*pi/180;

th2=3*pi/180;

a1=[1];

a2=[1];

a=[1];

for i=2:M

a1=[a1 exp(-1j*i*pi*sin(th1))];

a2=[a2 exp(-1j*i*pi*sin(th2))];

end

A=[a1' a2'];

Rss=[1 0;0 1];

Rrr=ARssA'+sig2*eye(6);

for k=1:180;

th(k)=-pi/6+pi*k/(3*180);

clear a

a=[1];

for jj=2:M

  a = [a exp(-1j*jj*pi*sin(th(k)))];

end

P(k)=real(1./(conj(a)*inv(Rrr)*a.'));

end

figure;

plot(th*180/pi,10*log10(P/max(P)),'k')

grid on

xlabel('Angle')

ylabel('|P(\theta)| (dB)')

How is angle of arrival estimated?

The estimation of angle in linear array antenna is based on the relative phase variation between the target location and antenna elements. This phase variation in spatial dimension between antenna elements will provide you with the location of the target (in terms of angle).

Thus, the basic approach to obtain an angle through array antenna is to perform Fast Fourier Transform (FFT) over the antenna elements. This conventional method with provide you with the angle of target but with very coarse angular resolution (around 2/N radians) (angular resolution using FFT over linear array)

For obtaining the angle with higher resolution, we can use CAPON beamformer or Super resolution algorithms such as MUSIC, ESPIRIT etc.

In order to answer your question,

  Algorithms such as Capon and Bartlett use the position of the source
  for calculate the auto- correlation matrix which is further used to
  calculate pseudo spectrum

This is not entirely true, as the auto correlation matrix in spatial dimension is the input to the CAPON beamformer and till that point the algorithm does not know the location of the target.

The angle of arrival of the target is determined by taking possible set of direction of arrivals and then suppress the other direction which is not under consideration. One has to scan through all possible direction of arrivals (using array manifold matrix) and then with the spectrum estimate the direction of arrival (normally choosing the highest peak).

For your kind reference:

Algorithm for CAPON is mentioned below:

Take auto correlation (A) (among antenna elements)
Take inverse of auto correlation matrix (Ainv)
Create array manifold matrix (M)
Obtain capon spectrum using = 1/ ((hermitian(M) * Ainv * M), where hermitian operation is the conjugate transpose.
Finding the peak spectrum location (provides angle of arrival)

  How is the angle of arrival being estimated if the source position is
  required to run the algorithm in the first place?

So we do not need to know the source dimension, but we have to find it among the possible source position/direction (obtaining the maximum likelihood of the actual direction with possible directions)


Brute forcing an aperiodic s-shaped signal to be periodic by mirroring it before performing FFT I have an aperiodic signal of an s-shape in the time-domain that I would like to filter and calculate the derivative of using the FFT. The jump at the end of the s-shape has, of course, a strong influence on the obtained amplitude spectrum. I am trying to eliminate the influence of the jump without giving less weight to the beginning and the end of my measurement by applying a window function.
Hence, as a first shot in the dark, my idea was to "brute force" the signal to be periodic by mirroring the s-shape at the end of the measurement, yielding a "bell curve". This approach yielded significantly improved results in the filtered signal and achieved a good approximation for the derivative. I am quite glad that this approach works and it helped me make a lot of progress recently.
However, I am still puzzled about why this approach works in the first place. Shouldn't the mirroring operation on my input signal significantly influence the output? Is there a theoretical reason for this approach working that I have overlooked?
Thank you in advance for all your thoughts and explanations.

Edit: I included a simplified example of what I did to brute force the periodic signal. As can be observed, the orange curve (ideal measurement without measurement noise) has an s-shape and ends abruptly at a certain value. I mirrored the curve around its last value to force a periodic "bell curve".

Cascading the initial signal with it’s reverse does two notable things which facilitates the OP’s further processing:
First and most significantly it minimizes the end discontinuities which are visible by creating the equivalent periodic waveform. The Fourier Transform result for the finite duration waveform is equivalent to the transform for the same waveform periodically extended for all time (the non-zero Fourier results specifically). This visualization is useful to more intuitively understand the effect of the discontinuities in an equivalent continuous-time sampled context: very high frequency components are needed to transition in a very short period of time; these high frequency components alias back to lower frequencies for a sampled system and we can see the result as “spectral leakage”).
Second, as a note of interest, it creates a linear phase response in frequency. This is similar but not identical to the cascade of any FIR filter with its reverse filter: the resulting coefficients after the convolution of the two will be symmetric and the frequency response will be linear phase with a magnitude response as the cascade of the two filters (thus having a “squared magnitude” response).  Consider if t=0 were at the center of the resulting mirrored time domain signal: the resulting frequency domain signal as the Fourier Transform would be real (symmetric about the vertical axis at t=0 in the time domain, or f=0 in the frequency domain will always be real in the other domain). By then shifting it to a later time as a causal signal, instead of being real in frequency, the resulting response with the exact same magnitude has an additional linear phase added but is otherwise not changed (a linear phase in frequency is a fixed delay in time). An example of this that may be more intuitive is a similar example in the frequency domain and how the time domain would be real because of the symmetry: Consider $2\cos(\omega t)$ which is a real time domain signal and has two symmetric impulses in the frequency domain at $\pm \omega$, (and we see this with Euler’s formula with $2\cos(\omega t)= e^{j \omega t} + e^{-j \omega t}$). With just one of those frequency domain impulses at $+\omega$, the time domain waveform would be complex and given as $e^{j \omega t}$.
This second point is an outcome of the theorem that the Fourier Transform of a real valued sequence is conjugate symmetric (so equivalently the Fourier Transform of a conjugate symmetric sequence is real), and is just describing how the result becomes linear phase (real plus an additional time delay so then complex with linear phase). Most significantly the OP has reduced the aliasing effects from the otherwise very large discontinuity in time while the original time domain sequence is completely recoverable from the frequency domain, so less prone to distortion effects with any subsequent processing in frequency prior to transforming back to the time domain, and noted that this was done without modification to any of the original time domain samples (as would be done with windowing). This was feasible since the derivative of the original waveform was zero or close to zero at the edges, thus minimizing the discontinuity in the equivalent periodic waveform after adding the mirrored component.


MATLAB: Compute pole-zero diagram of FIR using only impulse response Given a discrete time impulse response $h[n]$ of a system, is there a way to plot its poles and zeros in MATLAB? The input impulse response can be variable, so I can't compute its transfer function before hand. I have seen several options where given the $H(z)$ you can plot the Pole-Zero diagram but couldn't find any which computes it using just $h[n]$.

Since your response is finite, your system corresponds to a FIR filter, where the coefficients are given by the values of the impulse response. Hence, you can perform Z-Transform of the FIR and find its zeros. Finally, since you would just need to calculate the zeros of the nominator polynomial. See this MATLAB code, following your impulse response example:

h = [1 0 -2 0 0 3];
b = h; % the FIR coefficients of a filter

roots(b)

zplane(b,1);

Output:

ans =

  -1.6379
   1.2438 + 0.4498i
   1.2438 - 0.4498i
  -0.4248 + 0.9309i
  -0.4248 - 0.9309i

So, the zeros are directly given by the roots. For the poles, you have (N-1) poles at the origin.

Also, have a look at e.g. this answer.


FFT Pitch Detection methods: Autocorrelation or other? I've seen many posts here that state that an FFT is almost all you need to get the fundamental frequency from an input stream.  The problem is that, with sound, this doesn't work; we need to look for periodicity as well.  I know that autocorrelation is useful in this arena, but even then, I am unfamiliar with how one would extract the fundamental frequency from the data after autocorrelation.  Also, I am aware that there are alternative methods as well.

I am not interested in the easiest solution; I am interested in 1: the most efficient solution and 2: the most accurate solution, especially if these are mutually exclusive.

FFTs alone are lousy at pitch estimation for many really common sound sources, including male voices, and low piano or guitar notes.  Thus, the comments that an FFT is all you need are false and misleading.  (Although FFTs can be used as components of a more accurate composite pitch detection algorithm.)

When using autocorrelation, after you remove or discount subharmonic lag peaks, the pitch frequency may be around or close to the reciprocal of the time lag corresponding to the autocorrelation peak.

Pitch often refers to the human perception of such, and thus an estimate sometimes needs to include psycho-acoustics effects and illusions.  There are multiple methods for, and lots of book chapters and research papers on this problem of pitch detection.  How well each algorithm works might depend on your particular sound source and requirements.

Pitch detection/estimation methods include lag estimators, such as autocorrelation, weighted autocorrelation, AMDF and ASDF; and frequency domain analysis methods based on initial FFTs, such as cepstral/cepstrum methods, and the harmonic product spectrum algorithm.  Composite methods, and ones that use some decision analysis are less easy, but can be more accurate or robust for some sound sources.  Look for YIN, RAPT and YAAPT.  This list is by no means exhaustive.


Video stabilization using Lucas-Kanade method I got an assignment in a video processing course - to stabilize a video using the Lucas-Kanade method. Assuming the matlab code I wrote for performing LK on 2 images works (i.e, I can get the optical flow matrices u, v that allow me to warp one image towards another) - how can I use it to perform a video stabilization?

As far as I understand, the main idea is to calculate the OF for each 2 adjacent frames, and to warp all the images towards the first one. that means that I need to keep the total warp that I accumulate during the whole process- and practically: saving total warp matrices u_tot, v_tot and add to them the u, v matrices I got from 2 adjacent frames, to perform a warp of the latter frame according to u_tot, v_tot and so on (notice the matlab code below).

Am I going in the right direction? I'd like to get an advise... thanks

matlab code:

EDIT: changes the code + added the WarpImage imp. (23.4.17)

function StabilizedVid =LucasKanadeVideoStabilization(InputVid,WindowSize,MaxIter,NumLevels,NumOfFrames)
    % InputVid is a VideoFileReader object
    % WindowSize - window size for the LK phase
    % MaxIter - maximal number of iterations for the LK phase
    % NumLevels - number of levels in the pyramid for the LK phase
    % NumOfFrames - number of frames to process

    reset(InputVid);
    StabilizedVid = vision.VideoFileWriter('StabilizedVid.avi');
    cnt = 1;
    if (~isDone(InputVid))
        frame1 = im2double(rgb2gray(step(InputVid)));
        first_frame = true;
    else
        error('no frames in the video');
    end

    while ~isDone(InputVid)
        cnt = cnt+1;
        frame2 = im2double(rgb2gray(step(InputVid)));
        [frame1, frame2] = sizeCheck(frame1, frame2, NumLevels);
        if (first_frame)
            u_tot = zeros(size(frame1));
            v_tot = zeros(size(frame1));
            step(StabilizedVid,frame1);
            first_frame = false;
        end
        disp(['processing frames no. ', num2str(cnt-1), ', ', num2str(cnt)]);
        tic
        [u ,v] = LucasKanadeOpticalFlow(frame1,frame2,WindowSize,MaxIter,NumLevels);
        disp(['time: ', num2str(toc)]);
        u_tot = u_tot + u;
        v_tot = v_tot + v;
        frame2_warp = WarpImage(frame2, u_tot, v_tot);
        step(StabilizedVid,frame2_warp);
        if cnt >= NumOfFrames
            break
        end
        frame1 = frame2_warp;
    end
    release(StabilizedVid);
end

%%%%%%%%%%%%%%%%%%%%

function I_warp = WarpImage(I,u,v)
    % I - image to warp
    % u,v - the optical flow parameters gotten from the LK algorithm

    [x, y] = meshgrid(1:size(I,2),1:size(I,1));
    % interpulate
    I_warp = interp2(I, x+u, y+v, 'cubic');
    % in case of NaN values - put the original values of I instead
    I_warp(isnan(I_warp)) = I(isnan(I_warp));

end

For video stabilization, you need to estimate motion vector with respect to a reference frame $F_{reference}$ (which let's assume it to be the first frame). Then your motion model should be determined, for simplicity let's assume it to be Translational motion. Therefore, we need to estimate motion in $x$ direction and motion in $y$ direction for each frame. Now you have two options to estimate the change. If you are interested to estimate it through optical flow you can do either these:

Estimate regular conventional optical flow from frame to frame
Estimate Optical flow of first frame to the $i$ $th$ frame.

  I accumulate during the whole process- and practically...

If you've chose the first, then Yes, you need to accumulate the optical flow, but I think it would be simpler to go for the second method.

As a hint, I suggest track a token or compute optical flow of a small patch of the image.

Edit:

I feel the problem is with your Warp function, the code below is part of video stabilization assignment in my class, try something like this sample. The translational parameters are estimated using token based tracking.

function [ ou_I ] = frame_translation_func( in_I,Tx,Ty,ou_RoWs,ou_CuLoMns,ou_layer)
ou_I = uint8(zeros(ou_RoWs,ou_CuLoMns,ou_layer));
for i2=1:ou_RoWs
    for j2=1:ou_CuLoMns
        if (i2 - Tx<1) || (j2 - Ty<1) || (i2 - Tx>ou_RoWs) || (j2 - Ty>ou_CuLoMns)
            ou_I(i2,j2,:)=0;
        else
            ou_I(i2,j2,:) = in_I(i2 - Tx,j2 - Ty,:);
        end
    end
end
end


STFT - DFT size of the bins Having some trouble understanding this concept, really could do with some advice.

I want my STFT to have the following parameters: NFFT = 256 overlap/hop = 128

Now essentially, the algorithm will work as follows:

1) Split the signal into blocks of 256

This will result in around 72*265:

2) For each of these blocks, calculate against the Hanning Window

3) Create a "slider" that does an overlap of 128.. So, in theory it would
   therefore be: `size = 256 + 128`

Therefore, when I'm computing the DFT for each of the overlapped blocks, will my FFT remain at size 256 or will the size be 256 + 128 If this is the case, does each block in the resulting vector still have to be of size 256?

Thanks

EDIT:

This is now my result:

But, compare this to a spectrogram in matplotlib

Where am I actually going wrong? I cannot make any sense of this.

I've looked bat through all the data, the blocks are correctly overlapping, Hanning window is being applied correctly.

Could it be to do with the fact I'm using a 1D DFT? I.e.

std::vector<complex> FFT->transform(complex_vector[0...n], 256);

The overlap is the way that you walk on your signal to be analyzed, lets go to an simple example, imagine that your signal X has a length = 14, NFFT = 4 and overlap/hop = 2

X = 1 2 3 4 5 6 7 8 9 10 11 12 13 14

To know the number of blocks of analysis you will:

int(n_samples/Overlap) - 1

Then for this example:

(14/2)-1 = 6

lets get 4 points from X to pass to your FFT

At the first time you get from X

X2 = 1 2 3 4

from now you need overlap by 2

X2 = 3 4 5 6

Third time

X2 = 5 6 7 8

fourth time

X2 = 7 8 9 10

fifth time

X2 = 9 10 11 12

sixth time

X2 = 11 12 13 14


How to correct for delays in a PAM-X signal with orthogonal CDMA? I have a receiver which receives data bits at a rate of 5kbit/s. The number of data bits per burst is low (say, 16 bits). These databits are obtained after despreading the incoming signal using an orthogonal CDMA code (8-bit Walsh Haddamard code). This implies that the incoming chip rate is 5*8=40 kchips/s. The actual signal I get from the ADC is an upsampled I/Q version at complex baseband, sampled at a certain sampling rate (e.g. 5Msps).
I implement the system in gnuradio (but the question is general).

I figured that instead of explicitely performing the despreading via multiplication and summing up the values, it is easier to use the cdma_code as a filter. This can be interpreted as a form of matched filtering (I think). However, this simplified version misses an important point: Timing. In the worst case, the downsampling happens at the symbol or bit transitions and the result is garbage. Hence I use clock synchronization based on polyphase filtering (https://wiki.gnuradio.org/index.php/Polyphase_Clock_Sync).

However, I have a multi-stage downconversion system here: I am able time-synchronize the chip stream with the Polyphase Clock Sync block. However, when I apply despreading (by filtering the chip stream with the time-reverse chip sequence and decimating) I again have a time alignment issue: the 8 bit chip code has to be aligned with the input signal.

I tried using the polyphase clock sync block where the matched filter combines the matched filter and the CDMA code as suggested in Marcus Müllers answer below. However, I think such a system is unable to find the correct offset.

Consider this simple MATLAB experiment which generates a sample stream encoding $[1,-1,1,1]$. The receiver simulates the "Polyphase Clock Sync" block in gnuradio by finding the time shift that minimizes $(b\cdot db/dt)_{t=nT}$ where $b$ is a shifted version of the input signal:

fs = 400e3;
fchip = 40e3;
% test bits 1,-1,1,1  and add zero-padding
testdata = [ 0 0 0 0 0 1 -1 1 1 0 0 0 0 0 ]';
% 8-bit Walsh-Haddamard
H = [ 1, -1, 1, -1, -1, 1, -1, 1 ];
SF = size(H,2);
Nbits = size(testdata, 1);

% the test data at the chip rate
testdata_cr = upfirdn(testdata(:,1), ones(1, SF), SF);

% the modulated chip stream
chips = testdata_cr .* repmat(H(1,:)', Nbits, 1);

% upconvert chip stream to sample stream using rrcos (for simplicity)
hup2 = rcosdesign(0.8, 15, fs/fchip);
samples = upfirdn(chips, hup2, fs/fchip);

%% RECEIVER
% define the matched filter (=rrcos + CDMA code) at the receiver
mfilt = upfirdn(H(1,:), hup2, fs/fchip);

x = filter(mfilt, 1, samples);

% This implements a simplified version of the Polyphase Clock Sync block
delays = -fs/(fchip/SF)/2:1:fs/(fchip/SF)/2;
res = zeros(length(delays), 1);
B = zeros(length(x), length(delays));
for m=1:length(delays)
    B(:,m) = circshift(x, delays(m)); % works because of zero padding
    %plot(b);hold all;
    v = B(:,m) .* gradient(B(:,m));
    res(m) = sum(v(1:fs/fchip*SF:end)); % sum samples
end
% find the delays where b * db/dt is approximately zero
zeroxing = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);
idx = zeroxing(res);

%plot(delays, abs(res), '-', delays(idx), res(idx), 'o');
% plot the candidates
plot(B(1:fs/fchip*SF:end,idx))

Since this block uses $b$ as an objective function (using the delay index that minimizes the sum of the samples of $b\cdot db/dt$), there must be a global optimum which finds the correct sequence $[1,-1,1,1]$. However, there are multiple optimum delays (values close to zero):

Now I plot the resulting bit stream for all the delays where the objective function is close to zero and only one shift gives the correct $[1,-1,1,1]$ pattern:

I can reproduce this experiment with a gnuradio model as well (https://www.dropbox.com/s/48fy30vwkivahqz/test_polyphase_clock_recovery_cdma_direct.grc):

My question is: How to implement such a system and how would one do this in practice?

Couple of observations:

you can of course also correlate with the already pulse-shaped CDMA sequence as a whole – it's a long filter, then, but that (in AWGN) will essentially give you the basis of a maximum likelihood receiver.
Commutativity of convolution: Doesn't matter whether you pulse shape filter -> CDMA sequence filter or CDMA sequence filter -> pulse shape filter. The result is the same (and the same as if you first convolved CDMA sequence and pulse shape, and the signal then with that – see 1.)

In either case, you'd only do one synchronizer.


How to get started with fingerprint recognition? I have started with Matlab and Digital Image Processing. I am trying to do Fingerprint recognition. I found many research papers regarding this but I am unsure where to start. Concisely, what is the most basic approach/algorithm for doing fingerprint recognition ?

This  may help you.
Also, everything related to fingerprinting in here.
I think is a good documentation to start with and view some source code. Nevertheless, you should learn some other algorithms before you start coding.


What should be the singal duration in time doamin to correctly measure its frequency using FFT I am generating signal consisting of 4 frequencies 380,460,500,580 Hz using 20500 as Fs. I generate this signal for different durations of 0.004s,0.008s,0.02s,0.04s.

Fs=20500;
values=0:1/Fs:duration;
sig1=sin(2*pi*FREQ*values); % FREQ can be 380,460,500,580
sig2 , sig3, sig4
sigt=normr(sig1+sig2+sig3+sig4);

When I take FFT I get different results for different time durations.

   NFFT=2^nextpow2(length(sigt));
   Y=fft(sigt,NFFT)/length(sigt);
   f=Fs/2*linspace(0,1,NFFT/2+1);
   figure, plot(f,2*abs(Y(1:NFFT/2+1)),'r')

for 0.004s, I get one peak at 480.5 Hz
for 0.008s , I get one peak    at 480.5 and two lower side peak at
320.3 and 640.6.
for 0.02s ,I get    four peaks at 360.4,440.4,520.5,600.6
for 0.04s , I get four Peaks at 380.4, 460.4,500.5,580.6

All frequencies have time period less than any of above durations. Is there any specific criteria that signal must be in some length so that FFT can detect its frequency.
Is there any other method exist which can measure frequency for such short durations?
Please share any knowledge.

For separating sinusoids (resolution of multiple peaks with clear gaps between them in an FFT magnitude result plot), you normally need an FFT length approximately double the reciprocal of the difference in frequency between the pair that is closest together.  This is roughly double the length over which 2 sinusoids are orthogonal, and is required in order to see the gap between, and to have 2 clear peaks instead of just one fat one in an FFT result plot.  This length increases if a non-rectangular window is used.


Why do we deal with the eigenvectors of the autocorrelation instead of the data itself? How intuitively to understand why eigenvectors of the autocorrelation matrix are used, but eigenvectors of the matrix constructed from temporal samples have no sense and aren't used? For example, in detection of a harmonious signal in additive noise.

Some "gut-level" reasons why it is better to work with the autocorrelation matrix instead of a matrix with your observations:

If you want to take into account all your observations and you have a lot of data, you'll end up manipulating (inverting, multiplying) fairly large matrices. If you work with the autocorrelation matrix, you "summarize" your data once (in a fairly efficient step requiring just a FFT and an inverse FFT), and from then, you just manipulate your autocorrelation matrix of size $P \times P$ where $P$ is your model order (for example for AR modeling or sinusoidal modeling).
With some data it just doesn't work numerically to use the raw observations because you encounter situations in which you have to deal with matrices which are not guaranteed to be positive-definite.

For example, let us consider two approaches to AR model fitting.

Direct usage of the data matrix

The empirical quadratic reconstruction error on your data is:

$$\epsilon = \mathbf{x}^T\mathbf{x} + \mathbf{x}^T\Gamma \mathbf{a} + \mathbf{a}^T\Gamma^T\mathbf{x} + \mathbf{a}^T\Gamma^T \Gamma \mathbf{a} $$

where $\mathbf{a}$ is the vector of AR-coefficients, $\mathbf{x}$ is your observations vector, and $\Gamma$ the matrix with your delayed observations. You need to find the value of $\mathbf{a}$ that minimizes this. After derivation and a bit of shuffling, your solution looks like this:

$$\mathbf{a} = - (\Gamma^T \Gamma)^{-1} \Gamma^T \mathbf{x}$$

And you are screwed because you have absolutely no guarantee that $\Gamma^T \Gamma$ can be inverted. In the process, numerically speaking, you had to deal with fairly large matrix products if you have a long sequence of observations.

Random process view

if you adapt a "random process" angle to the problem, the quantity you have to minimize (the expected value of the error) is:

$$\epsilon = r_x(0) + 2 \mathbf{r} \mathbf{a} + \mathbf{a}^T \mathbf{R} \mathbf{a}$$

And you end up with the more palatable solution:

$$\mathbf{a} = - \mathbf{R}^{-1} \mathbf{r}$$

With a solid guarantee that this will be computable because $\mathbf{R}$ is positive definite!

It looks like your problem is that of sinusoidal modeling (rather than AR modeling). There's a lot of hand-waving here, but what I have said about AR modeling and the hurdles of using the raw data matrix ; also applies to sinusoidal modeling - with eigenvalue decomposition being the problematic operation instead of matrix inversion.