Image classification using SIFT features and SVM I am hoping someone can explain how to use the bag of words model to perform image classification using SIFT/SURF/ORB features and a support vector machine?

At the moment I can compute the SIFT feature vectors for an image, and have implemented a SVM, however am finding it hard to understand the literature on how use the bag of words model to 'vector quantize' the SIFT features and build histograms that give fixed size vectors, that can be used to train and test the SVM.

Any links to tutorials or literature on the topic is welcome, thanks

If you could implement an SVM, you can quantize the features. :)

Typically the features are quantized using k-means clustering.  First, you decide what your "vocabulary size" should be (say 200 "visual words"), and then you run k-means clustering for that number of clusters (200).  The SIFT descriptors are vectors of 128 elements, i. e. points in 128-dimensional space. So you can try to cluster them, like any other points. You extract SIFT descriptors from a large number of images, similar to those you wish classify using bag-of-features. (Ideally this should be a separate set of images, but in practice people often just get features from their training image set.) Then you run k-means clustering on this large set of SIFT descriptors to partition it into 200 (or whatever) clusters, i. e. to assign each descriptor to a cluster. k-means will give you 200 cluster centers, which you can use to assign any other SIFT descriptor to a particular cluster.

Then you take each SIFT descriptor in your image, and decide which of the 200 clusters it belongs to, by finding the center of the cluster closest to it.  Then you simply count how many features from each cluster you have.  Thus, for any image with any number of SIFT features you have a histogram of 200 bins.  That is your feature vector which you give to the SVM.  (Note, the term features is grossly overloaded).

As I recall, there was a lot of work done concerning how these histograms should be normalized.  I might be wrong, but I seem to recall a paper that claimed that a binary feature vector (i. e. 1 if at least 1 feature from this cluster is present, and 0 otherwise) worked better than a histogram.  You would have to check the literature for details, and the details are important.

Edit: The Computer Vision System Toolbox for MATLAB now provides the bag of features functionality.


Correction of a signal through a transmitter I am inquiring as to a practical way to solve a problem I have.

Basically, I need to transmit a signal, $x[n]$, through a seismic transmitter. (It will go through a D/A, etc). The transmitter that we have has a certain frequency characteristic, in that, it passes certain frequencies out better than others. (Nothing shocking, but for the sake of simplicity we can assume it is linear across frequency).

My signal $x[n]$ is somewhat wide-band. Hence, at the output of the transmitter, it might become distorted.

I am wondering, knowing my ideal signal $x[n]$ and knowing my transmitter frequency characteristics, is there a way I can 'massage' $x[n]$ before putting it through my transmitter, such that at the output I get the original $x[n]$ that I wanted?

Like I mentioned I am seeking guidance on practical and actionable techniques for something like this.

If the context helps this is for a seismic transponder app.

Thanks!

For the sake of not putting more stuff in the comments, I'll try describe how I would approach the problem, although I confess that I've never actually done this myself.

You mentioned yourself, that for the sake of simplicity, we can assume for now that the system is linear, so that all frequencies are distorted in a consistent way. Perhaps it's non-linear as a function of distance between transmitter and receiver, but let's assume for now that this distance is fixed. I also assume that time-invariance property holds.

Linear time-invariant (LTI) systems are nice in many ways, one of which is that they commute (note that this is generally not true for time-varying systems). In other words, if $L_1$ and $L_2$ are LTI, then $L_1(L_2(x[n])) = L_2(L_1(x[n]))$ This means that if you can work out how to correct the received signal for distortions caused by your channel after you've received the signal, the same correction can be applied as pre-distortion before you send it out. Note all the assumptions we made, so if these definitely do not hold, perhaps you should clarify and we can start relaxing them.

Now to the topic of how in the world do you compute the correction filter. There are two equally straightforward ways, depending on whether you need it to be and FIR filter or if you can allow it to be IIR.

The first thing you do is measure the impulse response of the system. This can be done in many ways. For seismic systems I've heard people use chirp deconvolutions and wavelet based methods (I don't know much about exact ones). I would start off by sending a chirp through your channel, recording it on the other end, taking FFT's of original and received signals, dividing recorded by original and taking the IFFT. I know I will probably get lynched for suggesting this, but this is a first step. This may actually work okay, so try it. You can then start experimenting with more sophisticated ways of obtaining the impulse response.

Once we've obtained this impulse response, call it $h[n]$, we need to invert it. From basic LTI system theory, we know that inverting $H(z)$ is done simply with $H^{-1}(z)=\dfrac{1}{H(z)}$. However, there's a catch. $H^{-1}(z)$ may not be stable if $H(z)$, which is a polynomial, has complex roots outside the unit circle on the complex plane. This can be remedied by converting $H(z)$ to a minimum-phase system. This can be done using real cepstrum, and if you have MATLAB with the signal processing toolbox, the rceps command will do that for you without any additional code (look up the documentation). Now we have an IIR all-pole filter that solves the problem.

If IIR is not a good option and you need an FIR filter, Linear Predictive Coding (LPC) can be used to approximate this filter with an FIR filter. Describing it is beyond what I'm going to describe here, but you're welcome to post more questions on these topics if you like. This is more of a system overview rather than the full solution, but working out the entire thing without knowing constraints and requirements of your project is impossible.


resolving clipping audio issues I have implemented a pre-emphasis filter by the following (pseudo) code:

a = 0.5;
s1[0] = s[0];

for (n = 1; n < N; n++) {
  s1[n] = (s[n-1] * a) + s[n]
}

The problem I am finding is that due to summing, there is clipping happening.  I thought the solution would be to find the difference between the max value in s1[n] and 1.0 and then subtract s1[n] by that amount...

However, that results in my audio turning into complete garbage, which I am confused as to why........  Subtracting a constant amount from an entire signal should simply lower it's amplitude by that amount, correct?

Subtracting a constant value from all of your samples will just cause all the clipping to happen at -1, and will give you some nasty distortion. The proper way to remove clipping is to normalize the samples.

Keep track of the max value of s1 in your loop, then in a second loop divide all of of your points by the max value, i.e.

a = 0.5;
s1[0] = s[0];
max_s1 = 0;
mean = 0;

for (n = 1; n < N; n++) {
    s1[n] = (s[n-1] * a) + s[n]
    if (abs(s1[n]) > max_s1) {
        max_s1 = abs(s1[n]);
    }
    mean += s1[n];
}
mean /= N;
max_s1 -= mean;
for (n = 0; n < N; n++) {
    s1[n] = (s1[n] - mean) / max_s1;
}

EDIT: For bonus points your can also subtract the mean from your filtered signal prior to doing the normalization. This will ensure you get maximum level without distortion from clipping.


NCO loses lock after a while for certain values of alpha EDIT: So the transfer function for my system is

$$ F_{nco}(z)=\beta E(z)/(z-1)$$
$$ z\Theta_{nco} (z)=\Theta_{nco} (z)+\alpha E(z)+\beta E(z)/(z-1)$$

So I would like to obtain $\Theta_{nco} (z) / \Theta(z)$ transfer function cause I think that would allow me to calculate whic range of $\alpha$ would be good for my PLL to lock. The thing is I don't really know how to obtain $ E(z)$ from $e[n]= \arg \big\{ x[n] \cdot (y[n])^* \big\} $

I'm trying to code my own NCO in order to use it for my own Costas loop (since I couldn't get this to work Implementing Costas loop for 4-PAM) and so far the NCO is working good except for it losing lock after a while. These are the update equations I'm using:

$$\theta_{nco}[n+1]= \theta_{nco}[n] + \alpha e[n] + f_{nco}[n]$$

$$f_{nco}[n+1] = f_{nco}[n] + \beta e[n]$$

So I generated this input with $f=50 Hz$ and $f_s = 1 kHz$ during $0.2 s$ (that is, 200 samples):

$$\begin{align}
x[n] &= e^{j 2\pi n f/f_s+\theta} \\
  &= \cos(2\pi n f/f_s+\theta)+j\sin(2\pi n f/f_s+\theta) \\
\end{align}$$

Then I'm running a for loop to update the NCO, which initially has 0 phase and frequency ($\theta_{nco}$ and $f_{nco}$). The parameters being $\alpha = 0.1$ and $\beta = 0.0025$.

The way I obtain the error is

$$e[n]= \arg \big\{ x[n] \cdot (y[n])^* \big\} $$

being $y[n]$ the output of the NCO:

$$y[n]=e^{j (2 \pi n f_{nco}[n] + \theta_{nco}[n])}$$

The thing is, it locks quite fast but it loses lock towards the few last samples (at around 0.2 s) (see attached picture). If I increase the number of samples, it loses lock even more.

I managed to fix it by changing alpha and beta to $\alpha = 0.05$ and $\beta = 0.000625$

So basically what I'd like to know is how I can know which are the best values for $\alpha$ and $\beta$ since I'm just randomly picking a low number for $\alpha$ and then doing $\beta=\alpha ^2 /4$. Thanks!

I will assume that your PLL is purely digital. Is that correct?

I assume your PLL consists of 3 blocks

The first block calculates the phase of your incoming signal (the signal you wanna lock to).

The second one will compare the phase of the signal to the phase of your "internal oscillator" yielding the error e[n]

The third one is your "controller" which I assume is some kind of PI controller with alpha being the proportionnal gain and beta being the integral gain.

Your phase detection algorithm that detects the phase of your input signal has some kind of low-pass effect. This creates some kind of delay between the output your PI controller and the measurement. This delay puts a limitation on the values of your alpha and beta parameters. I.e. if you increase the values too much, your PLL can be unstable and you will lose the lock.

If you change the phase of your incoming signal, how many samples does it take for the change to be effective? This will tell you approximately what your delay is. From this delay, you should be able to find alpha and beta gains that will give you a good gain and phase margin. You should try at least 6 dB and 45 degrees.

edit : How you obtain the value y*? Is is some kind of delay chain? Hilbert Filter? I assume there is some kind of delay involved...

edit 2 : Your controller seems ok, your gains are low. Your bandwidth is low and thus the closed-loop should be stable. I suspect that your error happens when the phase measured by algorithm transitions between pi and - pi. Picture this case

Phase NCO = (pi - 0.001)

Phase ref = (pi + 0.001)

The phase difference should be 0.002 radian right?

However, phase ref will not actually be (pi + 0.001) but rather (- pi + 0.001)

Yielding error = (-2*pi + 0.002). Which is not the real error, this has the effect of suddenly increasing the bandwidth of your controller, which is something you do not want to do. This is probably why it is stable at lower values of alpha. The glitch would be there but not as big.

I'm assuming your algorithm yields a phase between -pi and pi, if the algorithm yields a phase between 0 and 2*pi, the thinking remains the same.


Reliably matching a degraded image with the original out of dozens of candidates I have several collections contaning low-resolution captures of unevenly-lit presentation slides:

and thirty to sixty original high-resolution slide images:

How would you design a similarity measure between the low-resolution captures and the high-resolution slides that would nearly always give the highest score to the correct image pair (or, analogously, a distance measure that approaches zero only for the correct image pair)?

Experimental setup

So far I tried a machine-learning solution to get a better feeling for how difficult the problem is. I assigned each image pair a feature vector consisting of:

various grayscale histogram similarity / distance measures (correlation, intersection, $\chi^2$ distance, and Bhattacharyya distance) applied to both the entire images and the image quadrants,
Pearson's correlation coefficient between Haralick texture features,
Hamming distance between various image hashes (aHash, pHash, dHash), and
Levenshtein distance between the image OCRs normalized by the maximum OCR length.

For the training part of a dataset, I take each feature vector $\vec x$, let $\vec a\vec x + q = 1$ for matching images and $\vec a\vec x + q = -1$ for non-matching images. Then I use linear regression to obtain $\vec a$ and $q$. Finally, for the test part of a dataset, I use $\vec a\vec x+q$ as a scoring function to retrieve a sorted list of high-resolution slides for each low-resolution capture and look at the rank of the correct result.

Experimental results

With a small hand-annotated dataset, this is getting me the average rank of 6, which is better than random draws, but not nearly good enough for an application that would overlay the high-resolution slides over live capture; such an application will require an average rank that approaches 1. There is an opportunity for the application to be smart and only look at a small window of slides around the current slide, but it still needs to correctly guess the initial slide.

As you might expect, the feature vectors are also quite slow to compute (several seconds per a feature vector with a naive implementation on a higher-end quad-core laptop), which makes the solution unsiutable for a real-time application.

Your question has many solutions which can be good or bad, so I am giving the solution I came up with, hope it is a good one.
Since your images are only different in resolution and color I think the best way to compare them is to make their size the same and then to get rid of color and shading differences only compare their edges not the images themselves. This way, you can use comparison criteria that you've tried now independent of color and shading variations. I simply used sum of differences between two images in the code below and for the images you provided seems it works.

Through this code I try to give you an idea idea on how it is working:

I1=rgb2gray(imread('Low.png'));
I2=rgb2gray(imread('High1.png'));
I3=rgb2gray(imread('High2.png'));

[h,v]=size(I1);
I2=imresize(I2,[h,v]);
I3=imresize(I3,[h,v]);

d1=edge(I1,'log');
d2=edge(I2,'log');
d3=edge(I3,'log');

score1=sum(sum(abs(d1-d2)))
score2=sum(sum(abs(d1-d3)))

Some low pass filtering before calculating the scores might help.


Channel estimation using bit reversed order of FFT matrix Is it possible to use the bit-reverse order of FFT instead of FFT to estimate the channel?
Based on my reading, I found that bit reversed order of FFT matrix can diagonalize the channel, but unfortunately it can't estimate the channel based on pilot data. It means, that when having $y = h*x$, where $h$ is a channel, $*$ is convolution operation, and  $x$ is a signal. Therefore, in order to get the frequency signal based on bit-reversed order FFT matrix,  $Y = F^b y = F^b (h*x)$. Therefore, the signal $Y = HX$ where the $H$ and $X$ are the frequency-domain channel response and signal gotten using $F^b$, respectively.
Based on my analysis, I got that $H = Y/X$, but when estimating $H$ based on pilot data taken from $X$. For example, $X_p = X_{1:4:end}$, the channel response $H$ cannot be estimated using $H = Y/X_p$. Why ? and is it possible to estimate $H$ based on $F^b$ and $X_p$ ?

So, according to your text,
$$ H := F^bh,$$
i.e. what you call the bit-reversed order FFT of the channel impulse response.
Since "bit-reversed order" just implies you're permuting outputs of the FFT, $F^b$ is just a row permutation $\Pi_{BR}$ of the standard DFT matrix $D$.
Since the "proper" DFT actually leads to a diagonal $H$ (as a linear channel never "mixes" any other frequencies in), you already know exactly what to do to convert your $H$ to a diagonal matrix:
Reverse the reverse bit ordering:
$$\Pi^{-1}_{BR}H=\Pi^{-1}_{BR}\left(F^bh\right) = \Pi^{-1}_{BR}\left(\Pi_{BR}Dh\right) = Dh,$$
which, in best classical OFDM nature, diagonalizes the channel.
(Using the full rank of any permutation matrix, you can also show this is the only way to diagonalize the result of the bit-reverse ordered FFT; putting this into a ray of optimism: Those who set out to invent OFDM will in expectation invent OFDM!.)
Channel estimation relates to the way you want to correct the channel influence. Hence, once you diagonalized your system, it becomes "normal" OFDM, and all the normal OFDM methods applies.


Multiply FS Coefficients by 0 or 1 to get Low Pass Filter? I have written a bit of code to upsample and interpolate a sample waveform in MATLAB. It's at the point where I have taken a signal, upsampled it 2x and filled in the gaps with 0's, and then I found the fourier series coefficients. I know I need to choose which coefficients to multiply by 0 and which to multiply by 1 in order to get a low-passed signal which I will then perform an inverse Fourier Transform on, and this will give me back my interpolated upsampled signal.

First of all - why do I need to low pass this signal? I suppose I kind of added period information my bringing in those 0's and upsampling the signal, but I don't understand how to find a cutoff for this ideal lowpass I need, and how to translate that cutoff frequency into 0's and 1's scaling my FS coefficients. I heard I need to multiply every other sample by 0, but I don't know why this would be correct.

Here is my code:

function [a] = UpSampler(x)
N=length(x);
Nu = 2*N;
xt=zeros(1,Nu);

for n=1:Nu
    if(mod(n-1,2)==0)
        xt(n)=x(1+(n-1)/2);
    end
    if(mod(n-1,2)==1)
        xt(n)=0;
    end
end

a = getDTFS(xt);

%multiply a coefficients by 0 or 1
%Perform Inverse FT to get interpolated signal

end

This is generally better visualized graphically. When you upsample a signal by inserting zeros between samples you get replicas of the original.

Let's take a simple bandpass signal and plot it. All I'm doing is taking noise and lowpass filtering it.

y = filter(fir1(255, 0.2), 1, randn(1, 1024));
plot(20*log10(abs(fft(y))))

Now let's plot an upsampled version. Usually L is used to denote the integer amount of upsampling in some texts, let's go with L = 2 & 4.

subplot(211)
plot(20*log10(abs(fft(upsample(y, 2)))))
subplot(212)
plot(20*log10(abs(fft(upsample(y, 4)))))

Notice how we now have 2 & 4 replicas of the original spectrum. This is effectively what happens when you upsample a signal.

Now the problem is I have spectral content in places that it shouldn't be. This is where lowpass filtering is used. In the time domain I'm filling in the zeros with interpolated data and in the frequency domain I'm trying to reduce that unwanted spectral content to finally arrive at these plots. Notice how the new filter cutoff frequencies are at fc/L since I've effectively increased the sample rate by fs*L. The original bandpass signal went up to 0.2*fs.

subplot(211)
plot(20*log10(abs(fft(filter(fir1(255, 0.10), 1, upsample(y, 2))))))
subplot(212)
plot(20*log10(abs(fft(filter(fir1(255, 0.05), 1, upsample(y, 4))))))


Isn't Gardner's algorithm and Early-Late gate the same thing? I'm reading about some symbol synchronization methods and I find these two (ELG and Gardner) very similar.

Early-Late gate takes 2 samples, one at $t=T-\delta$ and the other at $t=T+\delta$, then from that you get the timing error; Gardner takes also 2 samples for $\delta=T/2$, one at $t=T-T/2$ and the other at $t=T+T/2$.

I can't really see the difference between these two, although there must surely be a few. Could someone explain that to me?

Thanks.

The Gardner Timing Error Detector is diagrammed in the graphic below, where two samples per symbol are used, and the error is determined using Prompt*(Late-Early), and when synchronized the center sample (Prompt) will be midway between two symbols. In contrast, an Early-Late approach uses (Late-Early), typically on a correlated symbol response, and when synchronized the Prompt sample which is between Early and Late will be centered on the correlation peak.

Note that the Gardner Timing Error Detector forms a Timing Error Discriminator by approximating the derivative of the waveform at the assumed zero crossing and weighing that by the actual sample value at the crossing location. (The difference equation $y[2n+1]-y[2n-1]$ is an approximation of the derivative).

In contrast, a Early-Late gate approach approximates the derivative at the peak correlation of the data symbol. When timing is synchronized in a Early-Late gate approach, the two Early and Late samples will be a fraction of the data symbol on each side of the peak correlation. They could be spaced 1/2 a symbol, making the spacing the same as the Gardner, but the key distinguishing feature is the lock point of the center sample.

The functional operation of the Gardner as forming a discriminator by multiplying the waveform with its derivative is demonstrated by the two graphics below.

The functional operation of an Early-Late approach is shown below for the case of timing being synchronized, too late and too early. A GPS receiver is a good example where this approach is often used, where the output of the correlator of the received sequence with the GPS acquisition code used for every symbol would have the triangular response shown below, which specifically is the magnitude of the correlator versus timing offset with the code. Three correlators would be used, with the Early and Late correlator spaced approximately one half of the correlation width as shown in the diagram. The Early and Late correlators would be used to keep the Prompt correlator centered on the peak, which would be used as the recovered symbol.

Notice the key difference between the two approaches. In the Garder Loop the derivative is fairly constant in vicinity of the lock point, and it is the waveform itself providing the discrimination curve (change of error magnitude versus change in timing offset) - while in contrast with the Early-Late approach it is the derivative itself that provides the discrimination curve.


Why do we need to construct Gaussian pyramid using SIFT detector I am learning about SIFT detection and descriptor. I am slightly unsure about why a Gaussian pyramid is built for the image.

I do understand that within each octave, we are applying the Difference of Gaussian filter at different scales to the image and finding for each pixel location, whether it is a local maxima. We do this also for the image at a smaller scale. So now we have marked out the pixels in two different scales of image that have local extrema values.
What do we do with the local extrema values in these 2 different scales ? Do we compute SIFT descriptors at the 2 different sizes of image ?

Actually, the purpose of all this is to approximate a Laplacian of Gaussian!
This computation is part of the corner detection of SIFT. You can find corners by examining extrema of the Laplacian of Gaussians (2nd order derivative). You use Gaussians for denoising, and a Laplacian to find inflection points. However, it is classical to not deal directly with computing this second order derivative, but to approximate it instead by difference of Gaussians for simplicity. This is sometimnes referred to as Gaussian blob detection.
You can find some illustration of this process here.
With this approach your difference of Gaussians (DoG) will have a fixed (small) size in pixels, so you're bounding the computational cost. But you still want to detect blobs of different sizes (because objects have varyng size,s, and the camera can be closer or further away from the scene), and you achieve this by "zooming out the image" with the various scales of the pyramid and keeping the DoG size fixed, instead of keeping the original image size and changing the size of your blob detector.


How to describe voltage fluctuations I am looking for a way to describe the properties of voltage fluctuations - is the RMS amplitude a suitable measure? Attached is a screenshot showing the random current fluctuations (in the case of strong noise this corresponds to a random walk)

EDIT:This random walk is the voltage fluctuations in a programmed neuron model. The "resting voltage" is $-65mV$ - this can also be seen as the "starting point" of the random walk. I would now like to be able to measure the "magnitude" of the fluctuations.
The data vector looks like this $v[i] = (-65,-65.12,-64.63.....$)
for i=1000 time steps, for example.
But if I calculate the RMS value, I get a value that is more or less the same size as the average (i.e. 64.3mV), which is nonsensical... I may also be measuring it wrong, but according to Wiki:
$$RMS= \sqrt{\frac{1}{n}\sum_{k=1}^n x_k^2}$$

Perhaps standard deviation is what you want. Standard deviation
$$\sigma= \sqrt{\frac{1}{n}\sum_{k=1}^n (x_k-x_{ave})^2}$$
(where $x_{ave}$ is the average of your sequence) quantifies the fluctuations of your sequence above and below the sequence's average value. The square of the standard deviation, $\sigma^2$, is the variance of your sequence which is a measure of the power of the fluctuations.


Group Delay of Mulitrate system how do I find the group delay of multirate system.(sample rate converter) , IIR non-linear phase filter
Say decimation by 2 (with lpf) = tf1(num,den)
then again decimation by 2 (with lpf) = tf2(num,den)
In total decimation by 4.
i can do the inidividual stages
gd1 = grpdelay(tf1,fvec,Fs)
gd2 = grpdelay(tf2,fvec,Fs/2)
but how do i find the combined one of the entire system
can two group delays be added ?
Edited :
Assuming the can be added, is the following scaling inside and outside the grpdelay function correct ?
gd1 = grpdelay(tf1,fvec,Fs)/Fs/2 % converting from samples to seconds
gd2 = grpdelay(tf2,fvec,Fs/2)/Fs/4
Then we can plot it
plot(fvec/,gd1+gd2);
xlabel('Freq MHz');
ylabel('Group Delay (s)');

If the group delay values you get are in seconds (rather than samples) then yes, you can essentially add the two group delays.
You can verify this by feeding an impulse padded with zeros into your system. You will get your filter response, which is then decimated and feed into your second stage. The second stage will also then apply it's filter and do it's decimation. Where the peak is in your output response is (relative to where the impulse is) will correspond to the group delay - this is of course assuming your filters have linear phase.
This is essentially determining an impulse response for your overall system, but your system is linear but it is not time invariant i.e. delaying your input by a sample will produce a different output not a delayed one, so it really isn't an "Impulse Response".


Offset effect in FFT magnitude i will use matlab code just for convenience, i get the same behavior using a C++ code.

I have two signals:

Fp=5;
A=0.00025;
Phi=0;
Off=0.5;
KADC_5VOLT=5.0/32768.0;

tc=250e-6;
tmax=30;
t=0:tc:tmax;

Buf1= int32 ((A*cos(2.0*pi*t*Fp+ Phi*pi/180.0)+Off)/(KADC_5VOLT));
Buf2= int32 ((A*cos(2.0*pi*t*Fp+ Phi*pi/180.0)+0)/(KADC_5VOLT));

Buf1=double(Buf1);
Buf2=double(Buf2);

performing a FFT on Buf1 i get amplitude 2.4246e-04 for the 5 Hz component ((this values seems to change randomly as a function of A, while A is small); while for Buf2 i get amplitude 2.6287e-04 for the 5 Hz component. If i set A > 0.0015, performing a FFT on Buf1 and Buf2 i get amplitude 0.0015 for the 5 Hz component;

if i use just:

Buf1= ((A*cos(2.0*pi*t*Fp+ Phi*pi/180.0)+Off)/(KADC_5VOLT));
Buf2= ((A*cos(2.0*pi*t*Fp+ Phi*pi/180.0)+0)/(KADC_5VOLT));

performing a FFT on Buf1 and Buf2 i get amplitude 2.5000e-04 for the 5 Hz component.

What is the reason of that? is seems that the problem is related to a quantization error, but why is the offset so important? Also: how can i get rid of it?

Since your offset is 3/5ths (plus some integer), the truncation fractions will be different, resulting in a different sum of quantization noise.  You can get rid of this difference by using an integer offset (including/after scaling) before truncation or rounding.  (e.g. not (2^N)/5)

Also, if your signal is not exactly periodic in your FFT aperture, then the result magnitude of very low frequencies will have a small dependence on phase, relative to the window (edge or center), and thus not exactly match the input amplitude.


Question about convolution of shifted function-simple I memorized the convolution formula, all good. I started doing examples and in one of them I had to calculate the convolution of a shifted function with another one.

I know the property $$H(x)=f(x)\ast g(x)\\H(x-k)=f(x)\ast g(x-k)=f(x-k)*g(x)$$

So I had no trouble finding the result but I had trouble writing the integral. Which of the following is correct?

$$\int_{-\infty}^\infty f(τ)g(t-k-τ)dτ\\\int_{-\infty}^\infty f(τ)g(t+k-τ)dτ$$

Considering an LTI continuous-time system with impulse response $h(t)$, its response to any valid input $x(t)$ can be computed with a convolution integral as:
$$y(t) = x(t) \star h(t) = \int_{-\infty}^{\infty} x(\tau) h(t-\tau) d\tau$$ and equivalently
$$y(t) = h(t) \star x(t) = \int_{-\infty}^{\infty} h(\tau) x(t-\tau) d\tau$$
due to the commutativity of the convolution operator.
Hence, shifting the output $y(t)$ by $k$ could be shown as
$$y(t-k) = \int_{-\infty}^{\infty} x(\tau) h(t-k-\tau) d\tau$$ or equivalently as
$$y(t-k) = \int_{-\infty}^{\infty} h(\tau) x(t-k-\tau) d\tau$$
Where, the first one equals $y(t-k) = x(t) \star h(t-k)$ and the second one equals $y(t-k) = h(t) \star x(t-k)$.


Use of the Dirac delta as a sampling operator As we are all taught in signal processing classes, the Dirac delta is not a normal function but a generalized function. It is sometimes defined as either the limit of some normal function whose support goes to zero as some parameter goes to zero, while keeping unit area. This apparently leads to the usual representation of the Dirac delta function as being defined and only having meaning under the so-called sifting integral,
$$
\intop_{-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-t\right)d\tau=x\left(t\right)
$$

Yet, pick up any book or paper on signal sampling and you will see the sampling process defined as a modulation process with the delta standing alone as a normal function, multiplying the analog signal to be sampled, like this:
$$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)
$$
Why is this legal? Why isn't the sampling process modeled as
$$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}\intop_{-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-nT\right)d\tau
$$
or
$$
x\left(nT\right)=\intop_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-nT\right)d\tau\text{?}
$$

Similarly, in the frequency domain, we know that a linear system's response to a sinusoidal input is a common measurement. Using the common definitions for $X$, $Y$, and $H$, we know that
$$
Y\left(f\right)=X\left(f\right)H\left(f\right).
$$
If
$$x\left(t\right)=\exp\left(j2\pi f_{0}t\right)$$
is a sinusoid, then
$$
X\left(f\right)=\delta\left(f-f_{0}\right)
$$
and
$$
Y\left(f\right)=H\left(f_{0}\right)
$$
where we have once again used the stand-alone delta as a sampling operator.

Can someone please clearly explain my confusion on this matter? Why do we seem to be so cavalier about how we use the Dirac delta when modeling sampling?

(While we're in the neighborhood, maybe someone could explain why we're even allowed to think about the Fourier transform of a sinusoid since it it is neither absolutely integrable nor square integrable.)

By the way, I find the discussion at Why is dirac delta used in continuous signal sampling? to be unsatisfying. Can someone do better?

The OP's "equation"
$$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)\tag{1}
$$
is incorrect; the sum on the right side of $(1)$ is the sampled signal $x_s(t)$, not $x(nT)$ at all.  That is,
$$x_s(t) = \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right).\tag{2}$$
This is not an ordinary function since it is a weighted sum of Dirac deltas, but if we pretend that it is indeed an ordinary function, then the  "value" of $x_s(t)$ at $t=mT$ is
$$x_s(mT) = \sum_{n=-\infty}^{\infty}x\left(nT\right)\delta\left(mT-nT\right) = x(mT)\delta(0),$$
while if $t$ is not an integer multiple of $T$, then the "value" of $x_s(t)$ is $0$. Note the quotation marks around the word value and avoid the usual nonsense of saying that $\delta(0) = \infty$ that is found in so many low-level DSP books.  $x_s(t)$ is not a function and it does not have any value at any real number $t$. What the sampled signal
$x_s(t)$ as defined in $(2)$ is is an impulse train; a collection of impulses spaced $T$ seconds apart with coefficient or amplitude of the impulse at time $t=mT$ being $x(mT)$, that is, the sampled value of $x(t)$ at time $mT$. As stated earlier, $x_s(t)$ is not an ordinary function of time; it is, as all the other answers have pointed out, a functional, that is, a function that operates on functions as a whole, not on individual specific values that the function being operated on takes on as ordinary functions do.   Now, it is a standard notion that
$$x(t)\delta(t-t_0) = x(t_0)\delta(t-t_0)$$ provided that $x(t)$ is continuous at $t=t_0$, and so if we re-write $(2)$  as
$$x_s(t) = \sum_{n=-\infty}^\infty x\left(nT\right)\delta\left(t-nT\right)
\tag{3}$$ we see that $(3)$ formally resembles the commonly used sum
$$x_m(t) = \sum_{n=-\infty}^{\infty}x\left(nT\right)g\left(t-nT\right)\tag{4}$$
which describes the modulation of a stream of pulses (of the form $g(t)$ spaced $T$ seconds apart) by the sample values of $x(t)$ spaced $T$ seconds apart. Small wonder then that textbook writers notice the analogy between $x_s(t)$ as given in $(3)$ and $x_m(t)$ in $(4)$ and impulsively describe the process that leads to $(2)$ or $(3)$ as a modulation process instead of a multiplication process. (At this point, it is worth reflecting that several modulation processes (e.g. AM-DSBSC which can be described as $x(t) \mapsto
x(t)\cos(2\pi f_ct)$) are indeed just multiplication of a carrier signal by a baseband signal).

So, why is this rigamarole useful to anyone?  Well, let's look at the Fourier transform of $x_s(t)$:

\begin{align}
 \int_{-\infty}^\infty x_s(t) \exp(-j2\pi ft) \,\mathrm dt
&= \int_{-\infty}^\infty \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right) \exp(-j2\pi ft) \,\mathrm dt\\
&= \sum_{n=-\infty}^{\infty}\int_{-\infty}^\infty x\left(t\right)\delta\left(t-nT\right) \exp(-j2\pi ft) \,\mathrm dt
&{\scriptstyle{\text{sifting property of Dirac delta}}}\\
&= \sum_{n=-\infty}^{\infty} x(nT)\exp(j2\pi f(nT)) \tag{5}
\end{align}
which is the DTFT of the discrete-time sequence $\big\{x(nT)\big\}_{n=-\infty}^\infty$ of sample values of $x(t)$ spaced $T$ seconds apart. More generally, if we have a (linear time-invariant BIBO) filter whose impulse response (dirty word!) is $g(t)$, then applying the input $x_s(t)$ to this filter, we get that the output is
\begin{align}
\int_{-\infty}^\infty \left( \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)\right) g(\tau-t)
\,\mathrm dt
&= \sum_{n=-\infty}^{\infty}\int_{-\infty}^\infty x\left(t\right)\delta\left(t-nT\right) g(\tau-t)
\,\mathrm dt\\
&= \sum_{n=-\infty}^{\infty} x(nT) g(\tau-nT)\\
&= x_m(\tau)& \text{as defined in }(4)
\end{align}
So, Equation $(2)$ does have some utility whereas the jury is still out on the OP's $(1)$.

The OP asks:

  Why isn't the sampling process modeled as
  $$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}\intop_{-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-nT\right)d\tau ??
$$

Well, quite apart from the royal confusion that the OP is creating by using $n$, a parameter on the left side of the OP's alleged definition, as an index of summation, note that the sifting integral property says that the value of that integral is $x(nT)$.  Consequently, the OP's purported definition reduces to
$$x(nT) = \sum_{n=-\infty}^\infty x(nT),$$
which is nonsensical, and matters don't improve very much if a different symbol is used for the parameter on the left side because all one gets is
$$x(mT) = \sum_{n=-\infty}^\infty x(nT),$$
which suggests that $x(nT) = 0$ for all integers $n$ that are not equal to the chosen $m$.  That's no way to run a railroad!

Turning to the frequency domain, the OP asserts that

  If
  $$x\left(t\right)=\exp\left(j2\pi f_{0}t\right)$$
  is a sinusoid, then
  $$
X\left(f\right)=\delta\left(f-f_{0}\right) \tag{4}
$$
  and
  $$
Y\left(f\right)=H\left(f_{0}\right).\tag{5}
$$

While $(4)$ is correct, $(5)$ is not. Actually,
$$Y(f) = H(f)X(f) = H(f)\delta(f-f_0) = H(f_0)\delta(f-f_0)\tag{6}$$ which does not equal
$H(f_0)$ as the OP claims.  Indeed, in $(5)$ the left side is a function of $f$ while the right side is a constant so that the standard interpretation of $(5)$ is that $Y(f)$ is a constant, which in turn implies that
$y(t)$ is an impulse or Dirac delta!  Since $H(f)$ is arbitrarily chosen, what $(5)$ seems to be implying is that the response of any filter to a complex sinusoid is an impulse, which is manifestly nonsense.  What is true is that if we assume the system with transfer function $H(f)$ (a.k.a frequency response) has input $\exp(j2\pi f_0 t)$, then the output is
$$y(t) = \mathcal F^{-1}(Y(f)) = \int_{\infty}^\infty H(f)\delta(f-f_0)\exp(j2\pi tf)\, \mathrm df = H(f_0)\exp(j2\pi f_0t)$$
which is the epitome of the definition of the frequency response function:

  If for each real number $f_0$, the complex sinusoidal input $\exp(j2\pi f_0t)$ at frequency $f_0$ produces complex sinusoidal output $H(f_0)\exp(j2\pi f_0t) = |H(f_0)|\exp\big(j(2\pi f_0t+\angle H(f_0))\big)$ (which is also at frequency $f_0$ but has amplitude
  $|H(f)|$ and phase $\angle H(f_0)$), then $H(f)$ is called the frequency response function or transfer function of the linear time-invariant system.

In short, the OP's whole question is based on massaging "alternative facts" (that do not withstand scrutiny) to arrive at incorrect conclusions.


Any estimation algorithms that allow a user to indicate that certain Solutions are impossible? Variations of the Kalman filter and other algorithms are used for navigation and target tracking. Often times certain solutions might as well be deemed impossible based on known topography in the area. Without accounting for the topography, the algorithms I know of might say that a submarine is most likely under the sea floor, or that an aircraft has tunneled into a mountain. Preferably, an algorithm would give the most likely solution after treating certain solutions as impossible.

Are there any estimation algorithms that do this sort of thing?
If there are such algorithms, what are they?

This series of papers tries to enforce such constraints.

The constraints are on the state of the Kalman filter:


Complexity: Number of MAC operations on a DSP Suppose a complex $(N \times N)$ matrix $\mathbf{T}$, multiplying a complex vector $\mathbf{x}$ with length $N$.

On a given application that I'm studying, which is irrelevant for this discussion, I can get the same results of $\mathbf{Tx}$ by the following operation:

\begin{align}
\mathbf{W_1} \Re\{ \mathbf{y} \} + \mathbf{W_2} \Im\{ \mathbf{y} \}
\end{align}
where are $\mathbf{W_1}$ and $\mathbf{W_2}$ are complex $(\frac{N}{2} \times \frac{N}{2})$ matrices, $\mathbf{y}$ is a complex vector of length $\frac{N}{2}$, $\Re$ denote the real part and $\Im$ the imaginary part.

I'm trying to asses which one requires more operations in a real implementation.

My analysis is that the first ($\mathbf{Tx}$) requires $N^2$ MACs, while the second requires

$ 2(\frac{N}{2})^2 = \frac{N^2}{2}$ MACs, plus $\frac{N}{2}$ sums.

I'm not completely sure about this, neither about how a DSP deals with complex numbers.
Am I making any mistake?

I'd appreciate your help.

It appears to me that your analysis is basically correct, except that in the second case there are $\frac{N}{2}$ sums, not $N$ sums.

I'm not sure what you are looking for when you ask "how a DSP deals with complex numbers", but generally they are represented in one of two ways: as separate real and imaginary values, or as a "complex" entity (structure/object/whatever- depends on the language) that contains real and imaginary members.


Smoothing 3D Data for the Second Derivative I am not a signal processing expert and my feeble attempts at solving this problem have come up short.

I have a C++ application which is being fed regularly-spaced (in time) 3D position samples.  The samples are generally accurate and, if plotted, appear to be quite smooth.  However, when the second derivative is computed (to look at accelerations), there are occasional/rare spikes in the data which are unreasonable.  (I know the general range of possible accelerations.  And by "rare", I mean that over about two ours of data, I generally see less than a dozen spikes which are out of the realm of possibility.)  I want to do as little to the data as possible, but smooth it enough to remove the spikes.

My first attempt was to use a Kalman filter in our math library.  This removed the spikes, but changed the data far more than what I would really want and seems heavy-handed for what I need.

There seem to be a lot of options out there for 2D data, but I don't have enough signal processing background to understand the trade-offs or necessarily what might translate well to 3D data.  Any education or suggestions would be appreciated!

I have always found that a "Median Filter", https://en.wikipedia.org/wiki/Median_filter   to be effective on streaming data.
A brief intuitive description is below; but you must take steps to insure you restrict the bandwidth to the real data source you expect.  While you might not know the sources/spectrum/timing of noise and spikes (in which case the median filter is actually pretty good) you hopefully know the signals and characteristics of the signals you are looking for.  For extensive analysis, still iffy IMHO, along the lines of PDF (probability distribution function) and spectral analysis: there are some IEEE articles; but my analysis of most of the articles were that they weren't really useful in the face of unknown interference.  Aside from accommodating only the real signal bandwidth I wouldn't try z-s domain filtering without a really good reason.

------------------------ Median filter abbreviated description
Basically you pass the data into a shift register and pick the middle value.  Examples always show a sort process; but since you lose only one value and gain only one value during a shift, you can bypass doing full sorts each step
Having said that: if the spikes are in the input data then that is where this should be applied; but if the spikes are real "jerks" (i.e. occurring in the second derivative) then that is where it should be applied.
Since the derivative process amplifies higher order noises one has to careful.  If I was doing what you describe I would impose the median filter on the source data, first differences, and then on second differences (acceleration).  And then examine the results; look inside the algorithm and understand what data is tossed!
There is an alternate form: you examine a shift register of the data for min's and max's and just toss the extremes.
Both techniques incur some data loss that has to be taken into account.  You can interpolate, which is a little expensive, or just accommodate a lower sampling rate; say one in ten samples will always be tossed.


minimum-phase phase via Hilbert transform returned values Following my previous question: HRIR Minimum phase I managed to compute the minimum-phase phase of a FIR filter (in my particular case, HRTF filters).
However I am not sure of the phase values returned by my function.
The Python function is the following. It expects an arbitrary HRIR as input. I am using a 44100 long FFT to have the freqency bins perfectly aligned at integer frequencies at 44100Hz sampling rate (so for example a sinusoid at 1000Hz will have a magnitude of exactly 1 with no spectral leakage)
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
from scipy.fft import fft, ifft, fftfreq

def min_phase_coversion(HRIR):

    HRIR_fft = fft(HRIR,44100)

    xf = fftfreq(len(HRIR_fft), 1/44100)

    phase=np.angle(HRIR_fft)

    minimum_phase = np.imag(-hilbert(np.log(np.abs(HRIR_fft))))

    plt.figure()
    plt.plot(xf, phase)
    plt.plot(xf,minimum_phase)
    plt.plot()
    plt.grid()
    plt.show()

For the following plots I am using the KEMAR MIT dataset, in particular the plot have been computed for the HRIR at [azimuth=144°, elevation=30°, distance=1.4m]. The selected HRIR is the following (one per ear):

Once I perform my minimum-phase computation using the above Python function, this is what I get (phase for a single ear):

However, looking at the plot, I can't understand the values returned by
np.imag(-hilbert(np.log(np.abs(HRIR_fft))))

Is the returned phase unwrapped? Or am i doing something wrong? The original phase "wraps" as expected once it gets over +3.14 or under -3.14 (+/- pi), while the minimum one doesn't.

I computed and compared the minimum phase HRIR and the original one.
This is my final code:
def min_phase_coversion(HRIR):
    '''

    :param HRIR: the desired HRIR impulse response to convert into minimum phase
    :return: the minimum phase version of the original HRIR
    '''

    HRIR_fft = fft(HRIR,44100)

    #computing magnitude, tested with sinusoid, it works.
    xf = fftfreq(len(HRIR_fft), 1/44100)

    fft_magnitude = np.abs(HRIR_fft/len(HRIR_fft))

    phase=np.angle(HRIR_fft)

    minimum_phase = np.imag(-hilbert(np.log(np.abs(HRIR_fft))))

    plt.figure()
    plt.plot(xf, fft_magnitude) #computing and plotting magnitude
    plt.plot(xf, phase, label='original phase')
    plt.plot(xf,minimum_phase, label='minimum phase')
    plt.legend()
    plt.plot()
    plt.grid()
    plt.show()

    min_phase_HRIR = irfft(np.abs(HRIR_fft)*np.exp(1j*minimum_phase), 44100) # |H(w)|*e^(jPhi(w))

    HRIR_original = irfft(HRIR_fft, 44100)

    plt.figure()
    plt.plot(min_phase_HRIR[0:512], label="Min_phase", linewidth=0.5, marker='o', markersize=1) #ricostruisce bene
    plt.plot(HRIR_original[0:512], label="original phase", linewidth=0.5, marker='o', markersize=1) #ricostruisce bene
    plt.legend()
    plt.show()

    plt.figure()
    min_phase_fft = fft(min_phase_HRIR, 44100)
    min_phase_mangnitude =  np.abs(min_phase_fft/len(min_phase_fft))
    plt.plot(xf,fft_magnitude)
    plt.plot(xf,min_phase_mangnitude)
    plt.show()

    return min_phase_HRIR

The results I get seem reasonable looking at the plots:

and the magnitude is the same as I expected (the original one and the min phase one are perfectly overlapping):

However, I am not totally sure of the phase behaviour (as mentioned in my question) and its returned values.


Clarification on how bandwidth effects input signal I know that the system bandwidth is the maximum range of frequencies that a system can handle. Is the reason that the input signal gets distorted as the bandwidth decreases a result of fewer segments of the signal being passed through the system? Does this mean that if I have a bandwidth greater than $B_{h1}$ I would have an output with no distortion?

The system $H(f)$ appears to be a lowpass filter. As you state, the fact that the bandwidth decreases results in a more distorted output. However, I would not say that "fewer segments of the signal pass through the system".

The Fourier transform of the input in your image is a $\mathrm{sinc}()$ function. This function extends to $\pm \infty$, so to get absolutely no distortion you would need an infinite bandwidth. As this is not possible, then distortion will always appear. Depending on where we "cut" the sinc's "tails", the output will be more or less distorted. If we let "more frequencies" pass, then the output will look more like the input.

The square pulse has very sharp edges - it suddenly goes from $0$ to $1$. This speed translates to high frequency components. As $H(f)$ is a lowpass filter, high components will be filtered out. To what extent this will affect the signal will depend on the bandwidth. As it decreases, more high frequencies will be cut off, leading to a slower transition at the output, which leads to that smoothed version of the square pulse.


Importance of Discrete Cosine Transform over High Pass filter? I have recently studied few work where the low and high frequency components of an image is accessed using the discrete cosine transform (one such example is presented here).

Low and High frequency components of an image can be easily accessed by applying a high pass filter in transformed domain.

My question is that what extra importance does DCT have over the simple high pass filter in accessing high and low frequency components of an image. In other words instead of applying DCT can't we directly access these components using High Pass filter?

Yes, if you apply a filter, say a high Pass, you get the low frequency components of the image, but still in the image domain. You need a Transformation, like the DCT or the FFT to get to the frequency domain, to see the separate frequency components of the image.


Low pass filter with different non linear phases I want to design low pass filters with different non linear phases. Means the magnitude response of all low pass filters will be same but having different non linear phases.

Design your desired magnitude response
Make it minimum phase
Cascade with any type of allpass filter to create your desired phase response


How to convert A/D data to the needed representation I have an Analog-Digital Converter used for getting digital data from LoadCell scales.
In order to use this data in my formulas I need to convert that data to the millimeters.
But how do I do it?

Two methods:

(1) theoretical - get the specs for your load cell, analogue hardware, A-D, etc, and do the math

(2) practical - calibrate your system by applying a series of known accurate loads, logging the A-D output and then derive a calibration offset and gain (assuming everything is supposed to be linear)

Ideally you should do both, i.e. do (1) first and then do (2) to calibrate out any offset or gain errors in your system.


Is the Bilateral Filter a Solution of Some Variational Method? I've been watching these video lectures about variational methods in computer vision.
In one of such video lectures a it is explained for example that the Gaussian filtering is a solution of the diffuse equation (so Gaussian filtering is a specific class of diffuse filtering). I do wonder if the same applies to the bilateral filter, because if so it can means that all such filters can be derived by some functional to be minimized. This is because the diffusion equation is a solution of a variational problem like

$$
E(u) = \int_{\Omega} \left[(u-f)^2 +  \lambda\lVert \nabla u\rVert^2 \right]d\Omega
$$

Again, is the bileteral filter a similar case?

Yes indeed.
You may have a look on work (Paper) by Michael Elad which is called On the Origin of the Bilateral Filter and Ways to Improve It or Analysis of the Bilateral Filter.
You may have a look on the way he derives the Bilateral Filter from a Least Squares problem with weighted regularization on the Gradient of the image in the form of:
$$ \arg \min_{x} \frac{1}{2} {\left\| x - y \right\|}_{2}^{2} + \frac{\lambda}{2} \sum_{n = 1}^{N} { \left( x - {D}^{n} x \right) }^{T} W \left( y, n \right) \left( x - {D}^{n} x \right) $$
Where $ {D}^{n} $ is Image Gradient Operator (Raising by $ n $ simply shifts it) and $ W $ is the weights matrix of the Regularized Least Squares.


Deriving of magnitude, phase response and impulse response of transfer function I have a transfer function as
$$H(z) = \frac{-0.0625z^4 + 0.25z^3 + 0.625z^2 + 0.25z - 0.0625}{z^4}$$
I want to derive magnitude and phase response of this equation.

can someone help me from here I am stuck.
Thanks

HINT: The frequency response has the form
$$H(e^{j\omega})=a+be^{-j\omega}+ce^{-2j\omega}+be^{-3j\omega}+ae^{-4j\omega}\tag{1}$$
which can be rewritten as
$$H(e^{j\omega})=e^{-2j\omega}\big[ae^{2j\omega}+be^{j\omega}+c+be^{-j\omega}+ae^{-2j\omega}\big]\tag{2}$$
Now note that the term in brackets is purely real-valued. I trust that you can take it from here.


Can Convolution neutral network train(learn) separately (train different times)? I am new in Convolution neutral network(CNN). My question is, is there a way to let CNN train separately?  For example, at very beginning, CNN only need to learning hardwriting 0 and 1. After the training is good enough(0 and 1 only), now, I want CNN to learn hardwriting 2. So, since 0 and 1 had learned really good, is there a way to let CNN just learning hardwriting 2 only instead of learning everything (0 1 2) from scratch ??

You know, just like human brain, after you already know how to distinguish 0 and 1,  we just need to learn 2, then we do can distinguish 0,1,2. instead of clean our memory first,then learning 0, 1,2 at same time. Because that efficiency is so low. Please advise me. Thanks.

This is called curriculum learning , and not only this approach is possible but also it is shown to be quite effective.
A few years ago I also came up with similar idea. Basically the idea is, instead of flooding the network with different examples, we can categorize the training data from simple to hard and gradually feed data ordered based on difficulty.
After searching for a while, I encountered this paper by Yoshua Bengio, which i think is among the important ones on curriculum learning. Also check out
FlowNet 2.0 which is a network trained with this approach.


For an LTI system, why does the Fourier transform of the impulse response give the frequency response? I know that for a given system, the Fourier transform of its impulse response gives its frequency response.
I want to find where this property comes from, but haven't been able to find if it's a definition or if there's a mathematical proof available, for both continuous and discrete-time systems.

Let $h(t)$ denote the impulse response of an LTI system. Then, for any input $x(t)$,
the output is $$y(t) = \int_{-\infty}^\infty h(\tau)x(t-\tau)\,\mathrm d\tau.$$
In particular, the response to the input $x(t) = \exp(j2\pi ft)$ is
$$\begin{align}
y(t) &= \int_{-\infty}^\infty h(\tau)\exp(j2\pi f(t-\tau))\,\mathrm d\tau\\
&= \exp(j2\pi ft)\int_{-\infty}^\infty h(\tau)\exp(-j2\pi f\tau)\,\mathrm d\tau\\
&= H(f)\exp(j2\pi ft),\tag{1}
\end{align}$$
where $H(f)$ is the Fourier transform of the impulse response
$h(t)$.  In words, for an LTI system with impulse response $h(t)$, the
input $\exp(j2\pi ft)$ produces output $H(f)\exp(j2\pi ft)$.
This is precisely what we have as the definition of the
frequency response of an LTI system (call
this $FR(f)$ for now):

  for every frequency $f$, the response
  to $\exp(j2\pi ft)$ is $FR(f)\exp(j2\pi ft)$ which is just
  a (complex) constant times the input complex exponential.

But $(1)$ shows that $FR(f)$ is just $H(f)$, the Fourier
transform of the impulse response $h(t)$, which is what you
wanted to prove.


Ear Test using MATLAB I'm trying to test my ear frequency range by using MATLAB. The thing I could not understand is that I can hear the frequencies above 20 kHz. My sampling frequency is enough to satisfy the Nyquist criteria (80 kHz which is limit for MATLAB). So there is no possibility to have an aliasing. What can be the problem? Thank you for any help.
Here is my code:

K               = 1000; % kilo
sampling_freq   = 80*K; % 80KHz sampling for sin wave
sampling_period = sampling_freq^(-1);
f               = 20*K; %kilohertz
t               = 0:sampling_period:5; % 5 second sin signal
testing_sound   = sin(2*pi*f*t);
plot(t,testing_sound)
sound(testing_sound,sampling_freq);

First of all - don't use 80kHz (I don't know why you came up with this number as a limit for MATLAB, whereas it is 1e6 in fact) but something that is more reasonable and standardized, i.e. 96kHz (or 88.2kHz if you wish...). Another thing is that your sound card might not be able to serve such high sampling frequency. You should check that in playback options --> Advanced tab and set it correctly. Last thing - you should drop using obsolete sound and use audioplayer - unless you use some ooold MATLAB.


Initial conditions for systems described in state space - LTI or not? Suppose we have some system given by

$$\begin{aligned}
\dot{x}(t) &= Ax(t) +Bu(t) \\
y(t) &= Cx(t)+Du(t)
\end{aligned}$$

where $x(t)$ are the state variables, $y(t)$ is the output and $u(t)$ is the input. All matrices are constant. The same question applies to the discrete case

$$\begin{aligned}
x[n+1] &= Ax[n] +Bu[n] \\
y[n] &= Cx[n]+Du[n]
\end{aligned}$$

It is known that a system with non-zero initial conditions cannot be LTI. However, if $x(0)\neq0$, I don't see why the system above would not be LTI as it is expressed. As far as I know, if a system is expressed that way, it has to be linear and, as matrices do not depend on $t$, it should be time-invariant too.

So we have a system that has to be LTI as it is expressed in state space with constant matrices, but it can't be LTI because it has $x(0)\neq0$.

I can't see the mistake in the reasoning that leads me to this absurd contradiction. Can someone point it out?

If you look at Chapter 5 of:

https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-011-introduction-to-communication-control-and-signal-processing-spring-2010/readings/MIT6_011S10_notes.pdf

which is entitled "Properties of LTI State Space Models" , equation 5.33 doesn't seem to have a problem with initial conditions, or any other book (I stand corrected, there is one book) that I'm aware of.  Unless Oppenheim was touched with insanity, I'm inclined to accept his characterization that initial conditions don't disqualify a LTI system as "not linear" by his use of the term "zero input linear".

At the beginning of the notes, (and in Oppenheim and Shaefer 3rd edition) a LTI system is given as:

$$
y[n]=\sum_{k=-\infty}^{\infty} x[k] h[ n -k]
$$
 which doesn't require $h[n]$ to be causal or stable.  $x[n]$ doesn't have to satisfy $x[n]=0 \;\text{for} \; n <0$ .

There is emphasis  in the text that one needs to consider the entire history of $x[n]$, not just for $n \ge 0$.

let $$x[n] = \hat{x}[n]+\tilde{x}[n] $$
  where
$$ \hat{x}[n] = \begin{cases} x[n]\; \text{for}\; n<0 \; \text{and}\\
0 \; \text{for}\; n\ge 0 \end{cases}$$

and
$$ \tilde{x}[n] = \begin{cases}0\; \text{for} \; n<0 \; \text{and}\; \\ x[n] \; \text{for}\; n\ge 0 \end{cases}$$
 by linearity.

$$
y[n]=\sum_{k=-\infty}^{-1} \hat{x}[k] h[ n -k]  + \sum_{k=0}^{\infty} \tilde{x}[k] h[ n -k]
$$
if $y[n]$ is causal,
$$
y[n]=\underbrace{\sum_{k=-\infty}^{-1} \hat{x}[k] h[ n -k]}_{\text{zero input linear}}  + \underbrace{\sum_{k=0}^{n} \tilde{x}[k] h[ n -k]}_{\text{zero state linear}}, \quad n \ge 0
$$

The essential point is that initial conditions account for prior input.  Where $n=0$ is referenced for $x[n]$ is arbitrary, which is another manifestation of time invariance.  Initial conditions are not some arbitrary value vexing the system.  if $x[n] = 0 $ for $n <0 $ the initial conditions are zero.

Let's try something else. Let $z[n]=\tilde{x}[n+1]$ (advance by one sample) and with $\tilde{x}[n]$, the system was LTI without controversy.
But now,
$$
y[n]=\underbrace{ z[-1] h[n]}_{\text{zero input linear}}  + \underbrace{\sum_{k=0}^{n} z[k] h[ n -k]}_{\text{zero state linear}}, \quad n \ge 0
$$
and now we have an initial condition. A forward shift of 1 sample would make an LTI system nonlinear?

The logical fallacy at the root of question is to use the definition of zero state linearity and applying it to the zero input case.


RMS level from impulse response I am trying to quantify the average differences between several impulse responses over the frequency range and over octave bands.
So, I am doing (in Matlab):
NFFT = 2^13;
x_fft(1:NFFT,:) = fft(IRs, NFFT);
spData = (abs(x_fft(1:NFFT/2, :))*2)/NFFT;
refPa = 2*1e-5;
R = nchoosek(1:size(spData,3), 2); % Calculate all possible combinations
aC = 0;
for idxNote = 1:size(R,1)
    aC = aC +1;
    diffsRows(aC ,:)=abs((20*log10(spData(:,R(idx ,1))/refPa)) - 20*log10(spData(:,R(idx ,2))/refPa));
end

semilogx(freqs,mean(diffsRows))

I can see in a paper that this has been done by filtering the impulse response with an octave band filter, then by calculating the "RMS levels in dB" for each band of the impulse responses (IRs) and subtracting the levels from the levels of a reference IR. I believe I am doing the same in my code (mine is over the whole frequency range).
Can be the RMS of an impulse response considered as its amplitude in the frequency domain?
Is this the correct way to go?
Thanks in advance for your help!

Please refer to Parseval's theorem.
If your impulse response is $h(n)$ and its DFT is $H(k)$, then
$$\sum_{n=0}^{N-1} |h(n)|^2 = \frac{1}{N} \sum_{k=0}^{N-1} |H(k)|^2 $$
So the square root of the sum of the amplitudes in the frequency domain, divided by the number of points in the DFT should give you the RMS. It is common practice to calculate RMS on a frequency band basis to get more meaningful results.


What is formula for quantization step size? I have read from different web links,different formulas of quantization step size as shown highlighted in attached snapshot,but what is actual/correct formula?

The actual quantization step size depends on the implementation of the quantizer, and would be by definition the input voltage change that would cause the digital output to increase by 1 bit.
You can find the specifics for the actual device you are using on the manufacturer datasheet, such as this example from ADI's AD9459 ADC, where the input range is from $\pm V_{ref}$ where $V_{ref}$ is an externally supplied reference voltage, and OR in the table is an out of range indicator:

The actual threshold down to 0.5 lsb levels is of no practical use since the self-noise within most ADCs typically exceed 1 lsb (otherwise the core could be used as a higher bit converter). Datasheets will also specific ENOB (Equivalent Number of Bits) which is dependent on sampling rate and input frequency. This captures all noise sources and provides a metric of the same number of bits the device provides compared to a "perfect" ADC. Note that this precision is maximized with an input signal that is typically 1 to 2 dB below actual full scale where clipping would occur (as we approach the actual full scale, spurious signals increase due to non-linearities which degrades overall dynamic range and therefore precision).
That said it is sufficient in most cases to simply use the total input range (up to the backed off full scale) and divide by $2^{b}$ where $b$ is the stated precision in ENOB.


How can I objectively measure train noise from a sound recording? I have made an audio clip of some local train noise, and run it through a spectral analyser:
https://youtu.be/ily4VEoEuSc
Here's a screen-capture if you can't view that:

As I watch the graph while listening, I can see three or four lines at the bottom, which are between 500Hz and 4000 Hz, that start and finish when my brain subjectively says the sound is annoying. About two thirds in there's a high pitched line (around 10K), that's much less annoying but also of interest.
I am trying to come up with some objective identification and measurement of this to help monitor the problem. I'm a software engineer with an engineering background, but I don't have much DSP knowledge, and that's what I think would help here. I'm happy using DSP software libraries, and some applied maths.
So my main two questions are:

What methods can I use to measure those lines (the two lines around 500Hz and 1100Hz look to be the loudest)? I've looked a little into 'tone detection', but although the lines are there, it's quite fuzzy. Potentially machine learning could identify the pattern, but so far as I know at best it would be 'yes/no' - much harder to have a continuous output. That's not to mention the training data needed. One final thought is a bandpass filter around those lines, but I think they may go up or down with the speed of the train.

How can normalise the amplitude? I don't have a recording of before when I thought the noise was fine. I'm wondering if there's some objective method to capture 'background sound energy' (not the right term), so I can put a meaningful number on the amplitude  of those lines. I can then survey locals to get an idea of what threshold people start to get annoyed by it.

I know that noise is a very subjective thing. However, I've been able to collect contacts of 100s of locals who agree this is the particular sound that's annoying. Trying to put a number on it may be futile, but I'd like to know why if so.
Background
I have a local train noise issue - the good news is the local transport authority are able to rectify it by lubricating the tracks. However - if the lubrication runs out, it's then quite an effort organising enough people to complain so that it gets refilled. I'm hoping by recording this data it'll be much easier to justify.

One of the radio stations around here used to have the slogan "It may not be your favorite song, but it has a lot of the same notes."  A DFT analysis can tell you what notes are present, and at what intensity.  With calibration, you can even measure the real world level.  This doesn't get anywhere near psycho-acoustics though.
What you could possibly use it for it to define an objective measureable standard of when the greasing needs to be done.  The loudness of the noise is going to be dependent upon where you measure it, so calibration and absolute values aren't that useful.  What you need, as you alluded to, is comparable recordings of the sound when the tracks are well greased vs when they get annoying.
With those in hand there are several techniques you can use to make quantitative comparisons.  The difficulty for you is that your intended audience will likely not have the technical understanding to comprehend any sophisiticated analysis so you are going to have to try to reduce down to a simple number or two, or a nice visual image that expresses the data succinctly.  I think you already have the latter with your spectogram.  For the former, a simple sum of the magnitudes of the bins in the frequency range of interest should be adequate.
Now, if you can track such a metric throughout a greasing cycle or two, then you can establish objective criteria (by association) of when the noise has reached the same level of annoyance.  Just because you are a concerned citizen, maybe you could email this value to the railroad company for every train passing and they will gladly monitor your numbers, thank you for your service, and send somebody out every time it reaches a threshold value.
Dream on?  Yep.
If there is a local noise ordinance, you can perhaps work with your local politicians  rather than the transport authority.


Noise Peaks Near Nyquist Frequency I am working on an audio upsampling/oversampling project where I am testing different interpolating schemes in hardware and comparing measurements to my expected results. In one of my interpolating strategies, linear (first-order-hold), I am measuring unexpected FFT peaks centered around the Nyquist frequency.

My setup is: between every two digital samples of recorded audio (44.1kHz), 15 intermediate samples are computed along a line (fixed point, in an FPGA). The 16x oversampled data is sent to a D/A (which zero-order-holds, obviously) and measured on a scope.

The results: if for example I am playing a tone signal of 2kHz, in FFT on oscilloscope I see main fundamental, small peak at (22.1-2)kHz, and a small peak at (22.1+2)kHz - unexpectedly! I also see as expected images centered around 44.1kHZ, 88.2kHZ, etc (these are also small, on account of linear interp./ZOH of DAC). I don't have any explicit analog processing.

Theories: I honeslty don't know...but possibilities are:

scaling I do in FPGA to convert from 16 bit input to 20 bits for DAC
using fixed point arithmetic in FPGA (I don't ever overflow or anything, but it rounds everything so it's not exact to be sure)
technically non-uniform sampling (I space the samples I created pretty evenly in time, but not perfectly...the last interpolated sample (i.e., #15 is held a little longer than the rest)

Any ideas on what gives?

You're doing interpolation "wrong", at least in the sense I assume you want it to happen.

So what you'd normally do is just take every input sample, and insert 15 zeros in between each of them. Spectrally, the resulting signal would be your input signal, plus 14 images of it, spaced 1/15 of the output bandwidth.

Now, your approach of sample and hold interpolation is the same signal, i.e. the zero padded one, unshifted, plus shifted by one sample, plus shifted by two..., plus shifted by 14 samples. So let $x$ be the input $y$ be the zero-padded and $w$ be your repetition-interpolated signal:

$$\begin{align*}
x[n]&:=\text{input}\\
y[n] &= \begin{cases}
x[\frac{n}{15}]& \text{for}\, n\, \mathrm{mod}\, 15 = 0\\
0& \text{else}
\end{cases}\\
w[n] &= x\left[\left\lfloor \frac{n}{15}\right\rfloor\right]\\
&= x[n] + x[n-1] + \dots + x[n-14]\\
&= \sum_{l=0}^{14}x[n-l]
\end{align*}$$

Hence, the spectrum of $w$ is just that of $y$, superimposed with 14 phase-rotated versions of itself.

$$\begin{align*}
W[z]
&= \sum_{l=0}^{14}\mathcal{Z}\left\{x[n-l]\right\}
\end{align*}$$

In fact, you've built a signal whose 15 polyphase decomposition components are identical!

Now, that will lead to interesting effects. In contrast, if you did the usual interpolation, i.e. took my $y[n]$ and low-pass filtering it to its original Nyquist bandwidth (i.e. $\frac{1}{15}$ of the output rate) it, none of the  polyphase components of the resulting signal wouldn't be identical to the previous one -- they'd incorporate a phase shift that ensures continuity.


What is the difference between Filters and Windowing? I'm finding some difficulties to understand the difference between Filters, and Windowing.

accoriding to my understanding both filters and window functions will make changes to the input signal parameters, according to wiki,

Filter:

a filter is a device or process that removes from a signal some unwanted component or feature

Windowing:

a window function (also known as an apodization function or tapering function[1]) is a mathematical function that is zero-valued outside of some chosen interval.

What I'm trying to understand is this,

Where do we use these windowing functions & Filters in practical world?
Do I need to use both for better signal quality, if so which one come first(Filters => Window or Window => Filters)?

Filters and windows complement each other, in a way. This is due to the convolutional lemma of the fourier transform.

Let me expand on that. A filter would be a convolution of a signal with another. In a technical sense, a window has a fixed (supposedly short) size and is seen as "cutting out" a short part of the signal. Often, a DFT would follow.

In the fourier domain, convolution and point-wise multiplication trade places, so a convolution in the time domain is a pw. multiplication in the fourier domain and vice-versa.

A filter scales the spectrum with a known frequency response. A window will filter the spectrum instead.

Filters and window functions often (but not always) have different design goals. For example, in a filter you would have rather large bandwidths. Windows percieved as useful have a small bandwidth (often, just a few samples). You can talk about both kinds in the same terms, you just have to remember that both are used in different contexts.


Construct IFFT from only the major frequencies of an input signal? I have data which when plotted looks like the green line in the below image. The data is from the walking-motion of a reinforcement-learning-model in a simulation and therefore the values do not really correspond to time. But the  dataset is 10000 steps long. I used the steps as the time variable. This is just a part of the plot.

To me the signal looks like a periodic signal with a lot of noise. Therefore my idea was to use Fourier Transform to identify a few underlying main frequencies and then simplify the signal. Thus cancelling out the noise.

I create the FFT like this:

s = ppo2_df['obs_4'].values

fft = np.fft.fft(s)
T = 1           # sampling interval
N = s.size

# 1/T = frequency
f = np.linspace(0, 1 / T, N)

plt.figure(figsize=(20, 5))
plt.ylabel("Amplitude")
plt.xlabel("Frequency [1/Step]")
# plot only half the spectrum
barlist = plt.bar(f[:N // 2], np.abs(fft)[:N // 2], align='edge', width=0.001)

max_inds = np.argsort(np.abs(fft)[:N // 2])[-5:]
for i in max_inds:
    barlist[i].set_color('r')

plt.show()

Next I try to reconstruct the signal from the FFT using IFFT using the 5 strongest frequencies discovered. My issue is that the eventual plot does not fluctuate around 0 like the original plot but around 0.4 instead. I don't know the mistake I must have made:

# filter n max amplitude frequencies
max_freq = f[max_inds]
max_fft = np.zeros(N)
max_fft[max_inds] = fft[max_inds]

x = np.fft.ifft(max_fft)
fig = px.line(y=np.abs(x[:1000]))
fig.show()

Can you tell me what I did wrong? (e.g. I am unsure if i chose a sensible sampling interval T = 1)

EDIT:

The data can be copied from here:

https://justpaste.it/3st5t

Without having the original data it is difficult to confirm, but the DC offset is due to bin 0 in the DFT. Bin 0 corresponds to the DC term or average offset value for the sequence. The plot in the DFT result shows a large value for bin 0 while the time domain waveform appears to be closer to 0 average. I would need to see the original data to determine how the DFT bin 0 occurred.

Update: Upon reviewing the original data from the OP, the average of this data is 0.4123 which is consistent with the result achieved.


Spectral Leakage in layman's terms I'm trying to understand the concept of spectral leakage for the DFT, without going deep to the mathematical intricacies (it's for practical purposes). I've read from the book "Introduction to Digital Signal Processing Using MATLAB with Application to Digital Communications" that:

  One reason for using DFT is to determine the spectrum of a
  discrete-time signal. If a frequency in the input signal does not
  coincide with one of the DFT bins, then the magnitude of the DFT of
  that particular signal will not be an impulse-like in shape. Instead,
  it will spread over the entire DFT bins. That is to say that the
  energy of the input discrete-time sequence in a particular frequency
  will leak into other neighboring bins. It smears the impulse-like
  frequency over the entire frequency bins. Therefore, the corresponding
  spectrum is only approximate. This is what is called the DFT leakage.

However it seems not clear to me. Are you able to explain more clearly the concept of spectral leakage for the DFT ? Maybe with an example if you can

Start with what we know:

The DFT works perfectly for a single cycle waveform input. For a single cycle waveform, all harmonics have an integer number of cycles (by the definition of “harmonic”).
“Spectral leakage” occurs when the input is not a single cycle waveform. As such, at least one frequency component is not represented by an integer number of cycles.

One way to view the DFT is as an algorithm that probes all possible sinusoids of the input waveform, determining frequency and phase of each. Fortunately, we have a limited number of sinusoids to probe for, knowing they must be harmonics.

One way to probe for a specific frequency, if you know the phase to expect, is to multiply by a sinusoid of the same frequency you’re looking for, at the same phase. If you take the average value of that point-by-point multiplication, it gives a value proportional to the amplitude of the component you’re looking for. Against any non-matching harmonic, the result averages to zero. Otherwise, it’s proportional to the amplitude of that frequency in the input. This is true with any number of harmonics present in the input.

Here's an example, probing for the first harmonic. There average of the results of the multiplication, averaged, yields exactly half the peak of the target:

And here is an example of the target not being the same frequency as the probe—the result wave has positive and negative points that average to zero, indicating no match:

The shortcoming for this scenario is that you need to know the phase of the harmonic component in advance. However, an alternative is to probe twice, once with a cosine wave and again with a sine wave. Because of the phase relationship they have, between the two results we can determine both amplitude and phase of the input harmonic under test. Please see my article A gentle introduction to the FFT for the finer details. In the DFT, the real part of the complex “bin” value represents its cosine result, and the imaginary part the sine result.

But, if a frequency component in the input is not an integer number of cycles, it isn’t really a sinusoid and isn’t a harmonic—easy to see if you try to repeat it end to end. And it won’t line up with any of the probe frequencies, and the math can’t have the property of either exactly matching or fully cancelling the individual probes. If the input is very close to a single cycles sinusoid, for instance, but having more than a single cycle, the probe-and-average with the expected first harmonic will give a value a bit less then a perfect match, the subsequent harmonic tests won’t full cancel either. That’s the “leakage”, if you want to view it as the missed energy spilling elsewhere.


How to expand on and deepen a basic skill set for image/video processing? I didn't take any class on image processing but I studied the subject based on lecture notes and important books for last three months. So at this point, I would like to sharp my image processing skills. Of course, I am practicing a lot of image processing work but the work seems to be very similar.

As a programmer, I learned a lot from top coder or similar where you can compete each other for better programming. Is there any web site I can compete for difficult image processing problems?

I hope you realize (from the form of your question, I think you do) that three months of reading up about it doesn't really make you an image processing expert.

I am familiar with topcoder, but while the difficulty of the problems and approaches that you needed to use there might be similar to those of image processing problems, to design a usable image processing application: do problem specific research, implement your findings and new ideas, test it to get reliable results, you need much more than a couple of hours (or even a week as for the topcoder marathon match).

If you want to just learn more about the basic image processing tools I suggest you:

take all your literature, random-open a page and study and implement the concept explained there
take a computer vision library such as OpenCV, or work in matlab, take some images, and try using and displaying various image processing methods on them (from simple morphology to more complicated things like Hough line detection). Try and predict results, graphically display them, and understand what you get.

If you are interested in specific problem in image processing (e.g. I am doing Content based image retrieval - CBIR right now -- and by that I mean, almost for the last 5 months), and want to get better in this specific topic, then my advice would be something like:

search google scholar, Mendeley article database, IEEEXplore for relevant works on the subject
download way too many papers, maybe 3-5 as many as you are prepared to read. Skim abstract, intruduction and conclusion sections and filter out the articles that sound relevant and interesting. Give priority to current (recent) articles, as well as articles with current references.
read all of that. Read the most important ones more than once. The one that will be the basis of your work is expected to have margins filled with your scribblings and to look at least a little bit scrappy by the time you finish your 5th read-through.
implement that. Go from simple, straightforward implementation to the optimized implementation.
before you start doing this, ensure that you have a dataset for testing. Testing is not relevant if it is not done on a large enough dataset. Some works mention the datasets they are using, or their datasets can be found on affiliated university websites.
there are public datasets (e.g. this one for object classification) that also organize annual challenges where you can test your new ideas against many peoples new ideas (but that's way bigger than topcoder :D)
if you get decent results, that's good. If you get better results than state-of-the art using some of your new ideas, double check them. Then triple check them. And then publish a fancy article ;)

I'm sure you can find a middle between just learning something about image processing just for fun and doing real time-consuming research from what I wrote... Actually, here's an idea: hang around here and try to understand and help to solve other people's problems! All of them need research-think-apply steps, they just vary with the depth of the steps they require ;) In any case, I hope it helps.


How to determine the symbol rate of a 2-FSK signal? Let's say that I have a 2-FSK signal recorded so that the SNR is around 20 dB and the noise seems to be more or less white. The distance between tones is greater than the minimum distance needed for non-coherent reception.

How would I go about determining its symbol rate?

I had some ideas, but it was a while since I actually had contact with FSK, so I'm probably forgetting something very obvious.

My initial idea was, since I have a recording, to use say spectrogram in Matlab and see how long each of the tones last. The downside here is that I need to be reasonably lucky with my data. I couldn't set up spectrogram parameters properly, so that end of each symbol is clearly visible on the time-frequency plots.

My next idea was to see if using analytical signal and getting the envelope would help. I don't think I obtained anything useful with that.

Later on, I tried to generate a signal of the same frequency as one of the symbols and then run a cross-correlation. I sort of expected the xcorr to be high when we're on the symbol of the same frequency as the one I'm using and not so high when we're not. What I got as the output is a sinusoid, but I couldn't figure out how to interpret it.

One approach could be to perform an FM-detection step (e.g. an atan2() operation followed by a first-order difference) to transform the waveform to measurements of the approximate received frequency versus time. Your FSK signal should then look like a binary-modulated baseband waveform. Then you can apply a nonlinearity to the signal to induce a discrete frequency component at the unknown symbol rate. This component can then be identified by spectral analysis of the signal after you apply the nonlinearity. Many different types of nonlinearities can work for this, including absolute value, squaring, and magnitude-squaring detectors.

If you're doing a literature search, the terms you might want to use would be blind symbol rate estimation. That implies the estimation of the symbol rate without foreknowledge of the signal's structure.


what should be the min number of samples and sampling frequency while playing sound As i am using open Al library for generating a tune. I was able to generate a tone for single frequency.

When I was trying to generate tone for multiple frequncies, it suffered from the problem of clipping and noise. So please sugest me how to rectify the problem.

Your sampling frequency should be at least two times higher than the highest frequency in your signal. (Read about Nyquist frequency).
This has nothing to do with clipping.
I can imagine that your clipping occurs because you are adding multiple signals without scaling (normalizing) the amplitudes.


Cooley-Tukey FFT: You don't have to zeropad to a power of 2? Someone wrote "The algorithm that Cooley and Tukey presented in their classic paper (Math. Comp. 19 (1965), 297-301. http://dx.doi.org/10.1090/S0025-5718-1965-0178586-1) can be applied to any composite length.  The performance advantages are greatest for highly composite lengths, of which powers-of-2 are one example, and lengths of powers-of-2 result in other advantages on binary computers, so it has become a common misconception that the algorithm is only applicable to signals whose length is a power of 2."

Does that mean that when you DO use the Cooley-Tukey FFT You don't have to zeropad to a power of 2?
Take for example an image of 1920x1080. So, if you want to use the Cooley-Tukey FFT, you don't need to zeropad that to 2048*2048?

Yes.  If you can factor length N, Cooley-Tukey will give you an DFT computation that uses less multiplies than an unfactored DFT.  A power of 2 factorization requires less lines of text or shorter equations to explain.  But any small prime factors allows a straight forward extension to the power-of-2 explanation or proof of the algorithm.

Many newer FFT libraries (iOS Accelerate for one example) allow efficient composite FFT lengths with factors of 2,3,5 and maybe more, thus requiring less padding to get to the next useful composite length.


Why is circular convolution used in DSP? Why not linear convolution? Why are we using circular convolution in DSP?
What's the main solid reason for the use of it in digital
processing?
Why does the concept of circular convolution come more often than
linear convolution?

Given a discrete time LTI system with impulse response $h[n]$, one can compute its response to any input $x[n]$ by a convolution sum: $$y[n] = x[n] \star h[n] =  \sum_{k=-\infty}^{\infty} {h[k]x[n-k]} \tag{1}$$

Without anything further stated, above definition is for the linear convolution (aperiodic convolution) between $h[n]$ and $x[n]$, which are aperiodic discrete-time sequences of possibly infinite length, unless otherwise stated.  This is different from a circular convolution which is between two periodic sequences of period $N$, and computed over a single period.

You can compute a linear convolution in time domain by Eq.1, or in frequency domain using the following DTFT (discrete-time Fourier transform) property:
$$y[n] = x[n]\star h[n] \implies  Y(e^{j\omega}) = X(e^{j\omega}) H(e^{j\omega}) \tag{2}$$

DTFT is naturally related with the linear convolution, as it deals with theoretically existing aperiodic sequences that may extend from $-\infty$ to $\infty$ reflected in its limits of the defining sum:
$$X(e^{j \omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} \tag{3}$$

When you want to make a calculation by hand, using symbolic expressions for signals, such as $x[n]=a^n u[n]$ and $h[n]=b^n u[n]$, you can calculate the results in time or frequency domains as outlined above.

Also, when you want to compute the same result by using a computer, you can use the time domain approach based on an LCCDE recursion (for IIR systems) or direct finite convolution sum (for FIR systems), BUT the frequency domain approach won't work with DTFT; as it's mainly a tool used to develop the mathematics of signals & systems theory, and it's not suitable for digital computer implementations, as its variable $\omega$ is a real continuous number.

What is used instead is the DFT (discrete Fourier transform) defined as $$ X[k] = X(e^{j\omega})|_{\omega = \frac {2\pi k}{N} } \tag{4}$$

where $k=0,1,...,N-1$ and $N$ is the length of the DFT, which's then called as an N-point DFT of signal $x[n]$.

Eq.4 implies that the DFT sequence $X[k]$ is obtained as the uniform samples of the DTFT $X(e^{j\omega})$, which is a periodic function, hence DFT $X[k]$ is also periodic but we only consider its first period from $k=0$ to $N-1$.

Since DFT sequences are inherently periodic, then their convolutions will also be periodic (circular). Therefore, whereas a linear convolution between aperiodic signals $x[n]$ and $y[n]$ is implied by the I-DTFT expression $$y[n] =\mathcal{I-DTFT} \{ X(e^{j\omega})H(e^{j\omega}) \}$$ instead a circular convolution between two periodic sequences is implied by the I-DFT expression $$\tilde{y}[n] = \mathcal{I-DFT} \{ X[k]H[k] \}$$

So, given that we want to compute a linear convolution between two aperiodic sequences $x[n]$ and $h[n]$ of lengths $L_x$ and $L_h$ respectively, using frequency domain by their $N$ point DFTs, $X[k]$ and $H[k]$, we actually have to compute a circular convolution between the periodic extensions of the signals $\tilde{x}[n]$ and $\tilde{h}[n]$ of periods $N$.

The key is in chosing a proper length $N$ of the DFTs, which must be long enough to avoid any time domain aliasing of the sequence $\tilde{y}[n]$, returned by the IDFT computation: $$ \tilde{y}[n] = \sum_{r=-\infty}^{\infty} y[n-rN] \tag{5}$$

where $y[n]$ is the result of the linear convolution that would be returned by the theoretical inverse DTFT and $\tilde{y}[n]$ is the periodic result of the periodic (circular) convolution implied by the inverse DFT.

Note that if any one of the signals are of infinite length, then it's NOT possible to compute their linear convolution using the DFT approach, as $N$ would go to infinity, practically impossible.
The implementation of a linear convolution via DFT then has the following steps:

Choose N according to the following criteria: $$N \geq L_x + L_h -1$$ which guarantees an alias-free reconstruction of the inverse signal $y[n]$ from its DFT $Y[k]$ of the computed circular convolution via $X[k]H[k]$.
Compute N-point DFTs $X[k]$ and $H[k]$ of $x[n]$ and $h[n]$.
Compute $Y[k] = X[k]H[k]$
Compute N-point inverse DFT of $Y[k]$ to produce the output $y[n]$


Digital implementation of type I compensation network I'm looking for advice or references on how to implement a "type I" compensation network of digitally.  An analog example is shown in the following diagram:

I am familiar with the basics of digital signal processing, and have read about digital integrators.  Would it be possible to use something like a first order IIR digital integrator, and just flip the sign of the output depending upon whether the input is above or below VREF? The goal is to drive a PWM output of a microprocessor, so there may be a better way to do this in the digital domain than attempting to copy the analog implementation.  Any advice would be appreciated.

The integrator is unstable, so it will probably just blow up. There's such a thing as a leaky integrator that does an OK job of getting pretty close. Have a look at the algorithm listed here on how you would go about implementing it using code. The math is pretty straightforward.

The differentiator uses b1=1, b2=-1. The integrator is just the inverse, but what you'll find is you need to modify (decrease) the filter coefficient by a subtle amount to make it 'leak'. For the integrator use something like a1=1, a2=-0.99. The closer a2 gets to 1.0 the greater the risk you run on it going unstable. Decreasing a2 to something like -0.95 will result in something that looks less and less like an integrator. Play around with it.


Tracking a Predefined Curve in a Series of Images I have a series of images. I want to be able to draw a curve on the image and track the curve in the series. For example, in the figure below, the left image illustrates the original image at time 0 and the one on the right is the curve I want to track in the next images in later times $t>0$. In fact, the edges that are marked as red change with time. In other words, the region, whose a portion of its boundary is shown in the figure, gets enlarged and shrunken in the subsequent images and I want to be able to track this.

Could someone kindly introduce me to one of the methods that is used to track a predefined curves according to the movement of the boundary which is the closest to it?

From my experience, I have successfully utilized Leo Grady's Random Walks method for this. The code is also available here. It works very well and can easily be made to run in real-time depending on the contour and image size. You could watch the video from my implementation here. Note that though, it might perfom differently than the original.


Best window before FFT ? (for a signal consisting of 2 tones used for phase measurement) I am doing phase measurments by transmitting and receiving tones (100 kHz and 16 kHz simultaneously). I am transmitting the tones  and receiving them, applying Blackman Harris window and doing FFT for phase measurements.
I feel the Blackman Harris window is not optimal for a signal that has only 2 tones but don't know which windowing function is better for phase measurement for such a signal.
Any Ideas?

Using an FFT to measure phase for just two tones results in a lot more processing that doing the following alternate approach that can be either streamed or processed in blocks. No windowing is needed:
Apply the received signal as the input to two multipliers. Apply a normalized local copy of one tone as the second input to one of the multipliers and a normalized local copy of the other tone to the other. The low pass filter of each output will be an estimate of the phase given by$A\cos(\theta)$.
For full 360 degree phase resolution, use a complex local tone and two complex multipliers with two outputs (I, Q) with the phase determined using atan2(Q, I).   (This is essentially the operation in two bins of the DFT).