Features of a non-periodic signal for comparison I am looking for comparing a signal (256 / 512 point discrete | non-periodic ) (obtained after extracting some features from the captured image) with a database of around 5000 signals (256 /512 point discrete | non-periodic)(constructed with the same features as the one in the image). The features are angles and distances from an image captured. The signal is constructed by taking angular value on the x-axis and the corresponding radial distance on the y-axis (amplitude). See the two images.

This image has 15 non-zero points and can be considered as one of the 5000 entries in the database.

Second image - signal acquired:

This image is that of the signal acquired and as can be seen it has 12 non-zero points in the signal. This is actually the corresponding signal of the first image but with some missing points and some variations in the amplitude.

The above is just a example of one of the entries of the database containing 15 non-zeroes and the acquired signal containing 12 non-zeroes. There can be many entries in the database containing the same number of non-zeroes, but the positioning of the non-zeroes will, of course, be different. The number of non-zeroes in a particular entry of the database varies from 5 to 50.

Now consider this as a signal to be compared.

There are two bottlenecks - Firstly, the acquired signal contains missing points from its corresponding signal in the database. The missing points can be very high sometimes. This can be viewed as the correlation between two signals can drop from 1 to 0.5 sometimes. Still I can apply correlation and look for the correct signal match.

But, here comes the second bottleneck and that is I do not want to correlate the acquired signal with each and every entry of the database (because it takes a huge amount of time). Instead, I want to shortlist some entries first based on some feature that can be calculated from the signal. This feature can be mean, variance, standard deviation, 3 dB bandwidth, etc. However, this feature must be something that does not vary much in the case that some points are missed as I stated earlier. I also thought of applying FFT on the signal and using some feature such as mean of the absolute value, phase information, etc.

Eventually, I have tried mean, variance, std, power bandwidth, mean of the absolute value of FFT. However, the change is quite drastic and I cannot shortlist to a small number (by small number I mean around 1000). At the same time, as I mentioned earlier, the correlation also drops quite drastically in some cases thus failing to find a match. Do you have any suggestions on how to go about for this problem of comparison when the acquired signal has missing points from the corresponding database signal? Any idea is appreciated! Thanks!

It is just an idea but how about a mean of the terms in a sliding window.

x = zeros(1,250);
c = randi(50,1,15);
x(c) = randi(500,1,15);
y = x;
d = randi(15,1,3);
y(c(d)) = 0;

figure;
stem(x,'LineWidth',1); hold on; stem(y,'LineWidth',1);
legend('Original Signal','Signal with Missing Samples');

N = 5;
L = numel(x)/N;
[m, m1] = deal(zeros(1,L));

for i = 1:L
    mx(1,i) = mean(x((i-1)*N+1:i*N));
    my(1,i) = mean(y((i-1)*N+1:i*N));
end

figure;
stem(mx,'LineWidth',1); hold on; stem(my,'LineWidth',1);
legend('Mean of x', 'Mean of y');

As you can see the plots are very similar now. You can also check their cross-correlation:

figure;
stem(xcorr(mx,my),'LineWidth',2);

These are of course just ideas to get you started. Now you can play with the window size and instead of doing $1:N$, $N+1:2N$ you can try $1:N$, $2:N+1$ , etc. and replace mean with other statistical parameters to see which one works the best for you.

Hope it helps!


why is it so hard to simulate up and down position on headphones? i'm wondering how come simulating front,back,left and right source position on headphones is easy but up and down is not ? (maybe even impossible ?)

is it due to the shape of the ear ? but should we not be able to modify the input signal according to the shape of the ear and then simulate ups and downs movements ? why binaural impulses don't simulate correctly ups and downs on headphones ?

thanks

Jeff

Front/back is actually very hard too, at least for stationary simulation without head tracking.

The main reason is simple: Left/Right is done by looking at the differences between the two ear signals, i.e. interaural level differences and interaural time differences. If a source is located to the left, the sound will arrive earlier on the left ear and will be louder (depending on frequency).

In the median plane (front, up, back down), the ear signals are almost identical. The only information available are pinna cues, shoulder reflection, interaural changes due to micro rotations of the head, small asymmetries, contextual information, etc. . That's a lot less robust so the localization acuity in general is a lot worse. It's also harder to model since it depends a lot on the individual human and one-size-fits-all doesn't work well.


extracting ultra accurate phase and frequency with FFT with real world data I was doing some experiments in Matlab, and write this code to experiment with it.

I want to know if it's possible to extract accurate within -+0.01° phase info from samples of an ADC, the matlab FFT function says that we can.

But it's weird somehow:   It would get me very accurate results, is it doable in a real-world scenario?

I assumed a 64 Sample per Cycle in my Wave under test.

Here is the matlab code

t=0:2*pi/64:2*pi-2*pi/64;
R=sin(t+pi*178.03/180);
FFT=fft(R);
angle(FFT(2))*180/pi

After this code I start to do some measurements on my real hardware, and I have increased my sample rate to 128 samples per cycle.

I have made some measurements with my system, which is a Cortex M4@200MHz device and an external ADC, I have a precision source injecting voltages and currents to my ADC which is connected to  CT and PT's, the first 4 channels are CT’s and the last 4 channels are PT's.

I have injected some sine waves to the system from 40Hz to 400Hz and in the 50Hz zone I have incremented the frequency by 0.1Hz, also I have tested the system with various phase angles with the channels, the sampling frequency of my ADC is 6400 sample/sec and my ADC is a true simultaneous one.

Here is the results, when I extract the phases respect to the V0 channel and in the 50Hz region, with different phases the extracted phase is within good ranges, but as the frequency changes the error would increase.

I have only 128 samples in a full 50Hz cycle and I cannot increase my sample rate higher due to my system's overhead. I have to extract all the phase and frequency of the waves from this 128 samples, roughly 5ms, aside of many other things to calculate.

Feel free to play with these data to see if you could suggest a good enough solution to extract phase and frequency info with under 0.1 degree for phases and 0.01Hz for frequency.

And all the extracted waves are ready to be injected into matlab.

Here you can get the wave files

https://ufile.io/5b7y5

You only have 64 samples per FFT. That means your FFT bins are spaced f_sample/64; if your sample rate is 128·50 Hz, then your frequency resolution, without interpolation, is just 100 Hz. To achieve a frequency resolution of $0.01\,\text{Hz}=\frac{50\,\text{Hz}}{50\cdot 100} = \frac{\frac{f_{sample}}{128}}{5000}= \frac{f_{sample}}{640\cdot5000}$, you'd need an FFT of length 3.2 millions. You can't do that on an M4, obviously, so: the FFT is not a tool applicable to your problem here.

Two approaches:

Either, you first reduce the sampling rate (including sufficient analog anti-aliasing filters) or you decimate your signal (including sufficient digital anti-aliasing filters) to a rate that describes a much smaller bandwidth of potential frequencies that you actually care about. Afterwards, do the FFT. Notice that this always always trades length of observation for precision: it's really that simple, you can't get 1/100 Hz of resolution without 100 s of observation. If things change within that period: bad luck.

Or, you don't do that and instead use parametric estimators for frequency and phase. I kind of feel I recommend that here. Parametric estimators can be as simple as zero-crossing detectors, or as complex as Eigenvalue-based estimators, and what you use depends on a lot of factors that you don't discuss – it's hard to guess what you need.

However, I like to point people at Peter K.'s A Review of the Frequency Estimation and Tracking Problems. Many of these algorithms will deliver phase info.

Regarding an accuracy requirement of 0.1° or 0.01 Hz in general: You can't have arbitrary precision with limited observation if you have any noise in your signal. Estimation theory doesn't allow for that. I think you'll be running into theoretical limits rather fast, depending on which algorithm you end up (and using the lowest-variance method might not always be the best choice!).


Getting error in PSNR calculation using Matlab I want to find out PSNR and SNR values of two images A and B of same dimension in Matlab.I used the following code

[peaksnr, snr] = psnr(A,B)

and getting an error

Undefined function 'psnr' for input arguments of type 'uint8'.

Then I converted both the images into double as follows

A = double(A);
B = double(B);

and again PSNR is calculated.Now I'm getting an error

Undefined function 'psnr' for input arguments of type 'double'.

Why this is happening?Please help me to troubleshoot the error.

Why this is happening?

The answer is quite simple and the error message exactly tells you what is wrong: Apparently you do not have the function psnr() available on your system.

This might have many reasons. Most probably:

You are missing a Matlab toolbox that provides this function
You are missing a Matlab script that provides this function

If you have the script, then Matlab does not know about it. Check your Matlab path settings


Calculate Frequency/Period from Output of a Maximum Entropy Spectrum MATLAB function When I plot said function output, "p," it looks like a typical spectrum with an easily identifiable peak at a given ix on the x axis. The help rubric from the function is

SYNTAX: p = mesa(x,m,nfreq);
For a vector x, this function calculates a maximum-entropy spectrum
of order m. The spectral estimate is returned in the vector p, which
has nfreq points linearly spaced in the Nyquist frequency interval 0-.5.
The psd is normalized such that the mean square value of x equals the
integral of p from -.5 to .5, so sum(x.^2)/N ~= sum(p)/nfreq.
Mesa is based on the Burg algorithm, as described in Numerical Recipes
and implemented in their memcof and evlmem subroutines.

Written by Eric Breitenberger        Version 5/24/95

My problem is that the x axis is unitless, i.e. if nfreq = 200 then the x axis is 1 to 200 long. If the peak is at x axis ix = 22, for example, how do I calculate the measured frequency/period of the input signal "x"?

To give a concrete example:

the length of the input vector "x" is 20 because I believe (or know due to a controlled input experiment) this length vector contains one full cycle (noisy sinusoidal) of interest

order m is set to 10, for example, although this is incidental to my question as I can see from experiment that changing the value of m only affects the sharpness of the peak

nfreq is set to 200

and with these values the peak occurs at x axis ix == 22.

From this, what is the mathematics to get the measured frequency/period from the output of the function?

Plot of the output

for a known input, e.g.

input_x = awgn( sinewave( 20 , 20 ) , 10 ) ;

From the code on p.$7$ of this document it can be seen that the $N$ (nfreq in the code) frequency points are chosen equidistantly between zero and Nyquist (i.e., $f_s/2$, with $f_s$ being the sampling frequency) using the command linspace(). So the individual frequencies are given by

$$f_k=\frac{kf_s}{2(N-1)},\qquad k=0,1,\ldots, N-1\tag{1}$$

[Note that the indices $k$ in $(1)$ start at zero.]


Particle filter Tracking I am tracking a car using particle filter by making a rectangle around it, and I am using the state vector $[x, y, u, v, a, h]$, where:

$x$, $y$ is the position of the body
$u$ and $v$ are velocity in $x$ and $y$ direction
$h = \mathit{height}$, and $a = \mathit{width/height}$ , height and width of the rectangle

My state equations for:

$$x = x' + uT$$
$$y=y'+vT$$

$u$ and $v$ remain same, and I am adding some random noise to it.

But I am not able to figure out, how should I change my $a$ and $h$ value, what should be the state equations for that. I am using histogram matching as the likelihood function.

The easiest way to do this is to keep track of the parameters a and h from detection. How did you pick out a car from the scene? a simplified method would be: background subtraction, followed by clustering, data association, find centroid of clusters => x,y. Also store the bounding box of the cluster, easiest way to keep an accurate a,h.
But I'm guessing you want to start with initial conditions a0,h0 and based on x,y,u,v update a,h. The problem is that a,h cannot be completely known using only x,y,u,v,T. You also need information about the imaging system or the environment.
Assuming we use a Cartesian coordinate system x=horizontal, and y=vertical
there are a few different ways to view this problem. If you are tracking cars from a bird's eye view (think google maps for instance)

then for all cases

a = a'
h = h'

realistically there can be slight variations. for instance if the car is facing this way __ vs this way | the bounding box would change, but given a set of x,y,u,v we have no way of knowing the cars orientation. Concretely, if the car is headed in the same direction (red and green) the bounding box shouldn't change. If the car changes heading (blue car) the bounding box changes drastically,

h = width' and width = h'

recall

a = w /h

so

h = a'*h'
a = width / h = h' / (a'*h')

unfortunately, we have no way of knowing with absolute certainty that the car changed directions. We can make up some kind though. Since cars are non holonomic vehicles ( a car can move forwards and backwards instantaneously but not left and right instantaneously) we can come up with some criteria to figure out a cars direction. For the red car we see |v| > |u| and |v'| > |u'|, or for the green car |u| > |v|  and |u'| > |v'| but for the blue car |v| > |u| and  |v| < |u| in this case we could guess that the car changed directions and the bounding box needs to be rotated using our second set of equations. However, if u~=v all bet's are off and we can't safely say which direction the car is pointing.
More complicated still is an image with perspective. For the following examples please note the origins in my picture are wrong, we should use the standard image convention where (0,0) is the upper left corner and (m,n) is the lower right corner where m=number of rows and n = number of columns. For this example I will use a image with a single vanishing point outside of the image.  We can split this into a few cases
CASE 1: car is moves horizontally along x

assuming the field of view is small we get

a = a'
h = h'

why? the height will only change as an object gets closer or further from us, since it is moving only along x, v=0 so our height shouldn't change. our width also shouldn't change (again the viewing angle must be small) because if we have a car traveling at velocity u, all points on the car move the x direction with the same u. More concretely, if delta_x = (x(T) - x(T-1)), delta_x should be the same for all points on the car. If the front tire moves delta_x (orange in pic), so should the back tire(light blue in pic), otherwise the car isn't moving as one solid body. IF both points moved the same delta_x, they should still be in the same position relative to one another (purple line in pic) so that means width doesn't change. we already said h doesn't change if the car moves only along x, and not y, so a shouldn't change either
CASE 2: car moves vertically, along y

this is a bit more difficult. As far as I know there is no way to determine this without performing some kind of offline calibration. There are a few issues, first is that size alone can't determine distance in an image. Think of an image with a person, and mountains in the background. The person looks much larger than the mountains, but we know this isn't true, however, a computer would have more difficulty figuring this out.
still using the single vanishing point outside of the image. We need to know something about the image, lets say we know the street is 6 meters  wide (call this Xstreet). Now we measure the street at various y locations. at the top of the image Xstreet0 = 100 pixels. At the bottom Xstreet1 = 180 pixels. Lets also say that Xstreet0 and Xstreet1 are 100 pixels apart vertically. Even though the width of the street appears to be changing we know in reality the street is really 6m at both places. So we come up with some kind of equation to mirror our perspective. Remember this equation

y = Mx + b

lets use that, so to find our slope (we aren't really interested in b), which models our perspective

M = rise / run = (y1,y0) / (x1 - x0) = 100 / 80 = 1.25

that means for every 1 pixel we move downwards in the y direction our x pixel gets 1.25 times bigger! With this knowledge we can now do some modeling

h ~= h' (this is really difficult to determine, it gets into perspective and turns very messy, since we said there is only 1 vanishing point outside of the image we can  say this doens't change much)
width = ((v*T)/M) + width'

looking closer we see that

v*T = delta_y
delta_y = M * delta_x ==> delta_x = delta_y / M ==> delta_x = (v*T)/M

so the term

(v*T)/M

gives us the change in X (delta_x), we simply add this to the old X to get our new X. but we aren't done yet a = width/height not simply width! so lets keep going

h = h' + delta_h
width = width' + delta_width
a = width / height = (width' + delta_width) / (h' + delta_h')

recall h ~= h' so delta_h = 0

a = (width' + delta_width) / h' = width'/ h' + delta_width / h' = a' + delta_width / h'

well we just calculated delta_width == delta_x, we know h' so we have our final equation

a = a' + (v*T)/(M *h')

where M was our slope (perspective compensation) we calculated offline.
CASE 3: a movement in x,y
is really just a movement in y, we have already said any change in x doesn't change the parameters at all. only a change in y

h = h'
a = a' + (v*T)/(M *h')

This problem gets a little more complicated when the vanishing line is actually in the scene. In this case, we have to keep track of if we are below or above the vanishing point and signs become very important.
This method also assumes that the objects we are tracking are in a constant z-plane. If our scene has a 3d axis, in our example, the street would be at one z value and cars travel in the x,y plane. We assume the cars don't levitate up and down, or our assumptions about y movements become invalid.
I mentioned it briefly, but this also won't work with wide FOV (field of view)cameras because our assumption that moving along x doesn't change width or height becomes invalid. If this doesn't make sense, just watch some fisheye video clips, like this video around 0 minutes 50 seconds
conclusion
Overall, my methods are highly approximate, and only work under certain conditions (assuming a camera doesn't move often, they are fairly reasonable restrictions). The absolute best way to keep an accurate a,h would be to save the bounding box from the original car classification steps. Hope this helps someone in the future


Parseval's Theorem for discrete series I need to use Parseval's Theorem to find
\begin{equation}
S\:=\:\sum _{n=-\infty }^{\infty }\:\left[\left(\frac{\sin\left(\frac{\pi }{4}n\right)}{2\pi n}\right)\left(\frac{\sin\left(\frac{\pi \:}{6}n\right)}{5\pi \:n}\right)\right]
\end{equation}
I used the fact that
\begin{equation}
x_1\left[n\right]\:=\:\frac{\sin\left(\frac{\pi }{4}n\right)}{2\pi n}\:=\:\frac{1}{8}\operatorname{sinc}\left(\frac{n}{4}\right)
\end{equation}
and
\begin{equation}
x_2\left[n\right]\:=\:\frac{\sin\left(\frac{\pi }{6}n\right)}{5\pi n}\:=\:\frac{1}{30}\operatorname{sinc}\left(\frac{n}{6}\right)
\end{equation}
Therefore, the Fourier Transform of each would be
\begin{equation}
X_1\left(e^{j\omega }\right)\:=\:1\:\:\:\:0<\omega <\frac{\pi }{4}
\end{equation}
\begin{equation}
X_2\left(e^{j\omega }\right)\:=\:1\:\:\:\:0<\omega <\frac{\pi }{6}
\end{equation}
so using Parseval's theorem
\begin{equation}
S\:=\frac{1}{2\pi }\int _0^{2\pi }X_1\left(e^{j\omega }\right)X_2\left(e^{j\omega }\right)d\omega \:=\:\frac{1}{2\pi }\int _0^{\frac{\pi }{6}\:}\left(1\right)\left(1\right)d\omega \:\:=\:\frac{1}{12}
\end{equation}
Is this correct? am kinda doubting I did a mistake in the theory somewhere. Thanks in advance.

You are right, but you have mistaken the DTFTs...

Considering the Parseval's relation for real discrete-time sequences to be:

$$ \sum_{n} x[n]y[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(e^{j\omega})Y(e^{j\omega}) d\omega $$

We shall define $x_1[n]$, $x_2[n]$ and find their corresponding DTFT as:

$$ x_1[n] = \frac{\sin(\frac{1}{4} \pi n)}{2\pi n} \Leftrightarrow X_1(e^{j \omega}) = \begin{cases} \frac{1}{2} &,& \text{ for } |w| < \frac{\pi}{4} \\ 0 &,& \text{otherwise} \end{cases} $$

$$ x_2[n] = \frac{\sin(\frac{1}{6} \pi n)}{5\pi n} \Leftrightarrow X_1(e^{j \omega}) = \begin{cases} \frac{1}{5} &,& \text{ for } |w| < \frac{\pi}{6} \\ 0 &,& \text{otherwise} \end{cases} $$

Where we have used the DTFT pair;
$$ x[n] = A\frac{\sin(\omega_c  n)}{\pi n} \Leftrightarrow X(e^{j \omega}) = \begin{cases} A &,& \text{ for } |w| < \omega_c \\ 0 &,& \text{otherwise} \end{cases} $$

Hence applying Parseval's theorem we see that:

$$ S = \sum_{n} \frac{\sin(\frac{1}{4} \pi n)}{2\pi n} \frac{\sin(\frac{1}{6} \pi n)}{5\pi n} = \frac{1}{2\pi} \int_{-\pi/6}^{\pi/6} \frac{1}{2}\frac{1}{5} d\omega = \frac{1}{2\pi} \frac{2\pi}{6} \frac{1}{10}  = \frac{1}{60} = 0.01667$$

To check this result you can use the following OCTAVE/matlab code for getting the sum in the time domain:

L = 128;    % select large enough L for better approximation
n = -L:L;
x =  sin(pi*n/4).*sin(pi*n/6)./(2*5*pi*pi.*n.*n);
x(L+1)= (1/4)*(1/6)/(2*5);
S = sum(x)
>> ans = 0.01667

which verifies that our calculation was correct.


Why images need to be padded before filtering in frequency domain In image processing books , we are told that Images need to be padded while doing filtering in the frequency domain. Why we need that zero padding ?

Multiplication in the frequency domain corresponds with circular convolution in the spatial domain.  This means that without padding the image properly, results from one side of the image will wrap around to the other side of the image.

You can think of 2D filtering as a sliding window that is centered over each pixel in the image and the center output pixel is a weighted sum of the pixels in the window.  With circular convolution, when the window hangs over the right edge of the image, it is actually wrapping back around to the left side of the image.  This means that output pixels on the right edge of the image will be affected by pixels on the left edge, which is almost never what is actually desired.

Zero-padding allows space for this wrap-around to occur without contaminating actual output pixels.


In OQAM modulation, why does the up-sampling process be seperated into two parts? The OQAM modulation is implemented by up-sampling the symbols by two and delaying the imaginary part of each sub-channel by half a symbol. The signals are up-sampled by M/2 and convolved with the impulse response of the pulse shaping filter u[n].Why should this process being separated into two parts?

In the end, it counts that the transmit signal of a multicarrier-OQAM system equals

$$x(t)=\sum_k\sum_m j^{k+m}d_{km} g(t-k\frac{T}{2})\exp(j2\pi\frac{mt}{T})$$

where $d_{km}$ is the real-valued data, $k$ is the symbol index and $m$ is the subcarrier index.

However you achieve this transmit signal structure, your implementation is fine. Though, it is mostly done as you describe for convenience reasons. For sure, you can perform the upsampling by $M$ points and delay one part of the QAM symbols by $M/2$ instead of first upsampling by two and then delaying one part by one sample. However, for implementation reasons, e.g. delaying by $M/2$ samples requires more memory than delaying by only one sample.

So, the answer is that it's not an algorithmic requirement, but more a convenient way of generating the OQAM multicarrier signal.


Do we really need two cameras for 3D vision? I heard that we can see the world 3D or we can detect distance because we have 2 eyes. but I think we can do it with only one eye! (test it!)

So why we can do it only with one eye? is it possible to have stereo vision for a robot only with one camera?

First of all, our brain does not only rely on our stereo visionary system to estimate the depth. There are many cues in a image scence for depth estimation, of which stereo, vision belongs to a sub-type called Binocular cues.

Technically there are many other methods of depth estimation, like Structure from motion, Perspective,etc. Just take a look at here and here.


FFT for real signals I am new to signal processing and especially to FFT. I need to know:

if it's possible to calculate FFT for example for vector $$V=\begin{bmatrix}1&2&3&4&5&6&7&8\end{bmatrix}$$ without calculate imaginary part
and will it be reversible ?

I suggest you do some more reading about the Fourier transform, but the short answer to your question is yes you can calculate the FFT for an arbitrary vector and discard the imaginary part, but it will no longer represent the Fourier transform of that signal. You can then perform the inverse FFT on that purely real spectrum, but it will not give you the original signal back (since it was not the actual spectrum of that signal).

If your starting signal happens to have even symmetry, then the spectrum will be purely real anyway, so you will be able to invert the real part of the FT and get your signal back.


OQPSK vs QPSK envelope query We use OQPSK because of envelope of signal but my question is

"In Phase Shift Keying(or QPSK) the amplitude of the transmitted signal is constant then how the envelope will change in OQPSK"??

Note that a QPSK signal using a band-limited pulse has an envelope that passes through zero every time there is a phase transition of $\pi$. So its envelope is not constant if all symbol transitions are allowed.

Offset QPSK (OQPSK) doesn't have phase transitions of $\pi$. By staggering the $I$ and $Q$ signals by half a symbol interval, the maximum phase transition is $\pi/2$, leading to an approximately constant envelope. Slight droops in the envelope (occurring at phase transitions of $\pi/2$) can be eliminated by hard-limiting.

Note that with an appropriate choice of the transmit pulse, phase transitions in OQPSK can be completely avoided. This results in a modulation scheme with an exactly constant envelope. One example of such a continuous-phase modulation is minimum-shift keying (MSK).

For more information on OQPSK take a look at this question and its answers. Also browse this site for other questions on OQPSK.


Producing Colored Noise from a given PSD Data I need to estimate (extract/produce) the colored noise of a given PSD Data. I have the following procedure to get the desired results.

First of all I calculate the filter coefficiencts using firls command in MATLAB by implementing Helbert Transform firls(filter_order,norm_freq,signal,'hilbert'), the problem I face is to choose the right order for the filter, I am at the moment working with hit and trial method and rougly a value of around 400 is selected which fits quite well with the PSD curve.
Secondly the white noise is passed to this filter and further amplified with a power amp. PSD of the extracted noise is again calculted and matched with the input PSD.

I want to know if this way is alright or not, moreover for passing white noise through the filter I am using Simulink, is there a way not to use time varibales and directly calculate the colored noise. In other word if I have to get rid of Simulink block. I tried to work with lsim and I can get to colored noise but it's amplitude is lesser and I am not sure how can I amplify the output signal so that it matches the input PSD.

just a short flow:

INPUT PSD --> Filter Coefficient Calculated --> White Noise Passed Through Filter --> Colored Noise at a lower amplitude extracted --> Amplified --> PSD of Colored Noise Calculated --> INPUT PSD matched with OUTPUT PSD

Estimating the necessary filter order
There are methods for estimating the needed filter order for low/band/high-pass filters, but the fact that you are talking about "colored noise" suggests that you are not going for relatively simple frequency profiles.  Thus, the answer is "no", I am not aware of an order estimator that will help you.  The best generic advice I can give you is to go as high as you can stand performance-wise and hopefully that will be enough.

Filtered/Colored output has lower amplitude than input
It sounds like you know the frequency profile you want, and are creating a filter that produces that profile.  If that is the case, fixing this should be easy.  Do the following:

Calculate the average of the amplitude of the frequency profile you are creating.
Divide all of your filter coefficients by that average.

That should make the overall gain of the filter "1" with regards to white noise input.


Detecting presence of music in ambient sound I am hoping to create an application that would listen to ambient sound and detect if music is being played. It is not important to identify the music being played; just detecting that some music is being played is enough.

I looked around for existing solutions but couldn't find any. Does anyone know algorithms that I can use to solve this problem? If source code is available, all the better. I found a couple of academic papers and implemented solutions suggested in them. But the results I obtained were not satisfactory.

PS: i) It would be a bonus if the algorithm is not computationally intensive; if algorithm is completely in time-domain that would be wonderful. ii) It is okay if the solution is not very accurate; occasional false-positives are okay.

This is a tough problem. The music identification service Shazam actually has a database of acoustic fingerprints for songs. Music is composed of many variable acoustic components and there are a slew of algorithms that may be used to estimate several components. One simple example would be the zero crossing rate, which is a time domain algorithm. I would refer to this link on the subject for a brief overview of pitch detection algorithms and go from there.


Differences between two closed-loop systems What could I say about stability of those systems by looking at the scheme?

System 1

System 2

$y_r(s)$ is input, $a(s)$, $b(s)$, $p(s)$, $q(s)$ and $r(s)$ are some polynomials and $y$ is output.

I calculated tranfer $\frac{y}{y_r}$ and it's the same for both systems.

$$\frac{y}{y_r} = \frac{r(s)b(s)}{a(s)p(s)+b(s)q(s)}$$

What makes the difference between stability and behavior of those systems?

Some progress:

I've edited the schemes. There is something, that obviosly differs in those systems. It's $e_{ss}$ - stady state error.

For the first system is:
$$e_{ss1}=y_r(s)\frac{r(s)}{p(s)}-y(s)\frac{q(s)}{p(s)}$$

For the second:
$$e_{ss2}=y_r(s)r(s)-y(s)q(s)$$

Does it tell me something about stability?

I guess this form is more usefull:

$$e_{ss1}=y_r(s)\frac{r(s)a(s)}{a(s)p(s)+b(s)q(s)}$$
$$e_{ss2}=y_r(s)\frac{r(s)a(s)p(s)}{a(s)p(s)+b(s)q(s)}$$

The second system's steady state error has higher order numinator than denominator. Does it mean, the signal's higher frequencies would be infinitely amplified?

There are more than just one way to approach the solution:

1) You can find the steady state error by finding the equation of the denominator of the transfer function.

In system 1:

$$D(s) = a(s)p(s) + b(s)q(s)$$

The order of the system is determined by this order of $D(s)$. Solving for $D(s)$ will give you the poles of the system. If the poles 'solutions of $D(s)$' all are negative. That means if all the values of the solution to $D(s)$ are in the left half plane $<0$. This means the system is stable.

2) The error in time domain is defined as "input - output".
That is:

$$e(t) = y_r(t) - y(t)$$

Now, to compute the steady state error, that is time at infinity 'after a long time'.
The equation will be: $$ e_{ss} = \lim_{t \to \infty}  y_r(t) - y(t) $$
That is simply taking the limit as t goes to infinity.

Now, using the final value theorem to compute the steady state error in the $s$ domain the equation will be:

$$E_{ss} = \lim_{ s \to 0} E$$

That is in your terms

$$E_{ss} = \lim_{s \to 0} s\cdot e_{ss1} \quad(1)$$

$$E_{ss} = \lim_{s \to 0} s\cdot e_{ss2} \quad(2)$$

Therefore; you can note that the result depends on the 'Type of the system and the Input'.
Type of the system: is how many poles at zero do you have?
That is writing:

$$a(s)p(s) + b(s)q(s) = s^n (s+a_0)(s\cdot a_1) \cdot ... \cdot(s\cdot a_{n-1})$$

the value n represents the order of the system.

Type zero means no poles to zero.
Type I means one pole at zero, etc.

Now if the input is a unit step that is $$y_r(s) = \dfrac{1}{s}$$
this will cancel out with the $s$ as specified in equations $(1)$ and $(2)$. The equation will be:

$$ e_{ss1} = \lim_{s \to 0} \dfrac{r(s)a(s)}{a(s)p(s)+b(s)q(s)} $$

Now, looking at the type of system:

If it is of type one or higher, it is clear now that the $e_{ss1}$ goes to
infinity yielding unstable system.
If type zero, it will have a certain value, yet still stable.

These are the techniques, and the answer highly depends on the type of the system or even the characteristic equation that is the equation of the denominator of the transfer function.


Camera Calibration using Single Input Image Given a camera, permanently located at the same position, angle and distance from what is being captured, and a single calibration image ( a row of vertical lines [ | | | | ] ), how would one go about using the image to apply correction (barrel distortion) to all subsequent images?

Numerous examples exist online for the camera calibration problem, but most attempt to achieve the correction through guess-work, or highly involved processes to correct real-world images. However, I cannot find anything on calibrating a camera whose position and subject is always at the same position, and where the expected positions of vertical lines are known prior to correction.

The closest I've come to my approach is the plumb-line calibration method described by brown (1971).

Are there any algorithms or pseudo code out there which describes an approach to calibrate using a single input image of vertical lines (where distortion would be clealry visible by the gap between each line)?

Close-Range Camera Calibration - D.C, Brown

Take a look at this paper here: "Straight lines have to be straight" by Faugeras et al.
http://link.springer.com/article/10.1007/PL00013269#page-1

It is straight forward to implement, but essentially removes distortions by observing the straightness of lines.

It should be straight forward to implement if you remove some automation (e.g. you click on the lines manually rather than relying on image processing).


Is there a relation between abs(fft(u.^2)) and abs((fft(u)).^2 Is there a relation between
(1) abs(fft(u.^2))
(2) abs((fft(u)).^2
where u is a general variable.
I have accidentaly found that if (1) has a power law +2, (2) has a power law +1. Can someone help me?
Thanks

Is there a relation?

Generally no. Squaring is a highly non-linear operation and completely changes the amplitude statistics and spectrum of the original signal. A simple example: if your input consists of a few sine waves of different frequency, the output will be all positive and the spectrum will contain lines are all possible combinations of sum and difference of the original frequencies. The bandwidth will also double which in the discrete domain will often result in aliasing.

has a power law +1

I have no idea what that means.


ROC vs stability in z domain? I have read in some dsp texts that when ROC includes unit circle,system is stable

But i am bit confused in difference between stability and marginal stability depending Upon ROC

Especially i am confused by part c of below dsp.se question

Z transform stability

Where in part C, ROC is 1/3 < |z| < 1 and analysis answer of part c is written as "non-causal & marginally stable" . I am confused by marginal stability. why it is marginally stable?If our ROC is changed to 1/3 < |z| < 1.1, will now it becomes stable or will it still remains marginally stable?

Marginally stable means that an otherwise stable system has one or more simple poles on the unit circle (in discrete time), or on the imaginary axis (in continuous time).

The consequence of that is that transients don't decay, but they also don't grow without bounds. Marginally stable systems are unstable in the bounded-input bounded-output (BIBO) sense.


How can I experimentally find the bandwidth of my PLL? I recently designed a phase-locked loop. I would like to quantify its bandwidth so that I can better understand it.
However, I made this bandwidth in an informal manner due to my lack of formal education on the topic. Is there an experimental way for me to find the bandwidth of this system?
My PLL uses a dual-channel lock-in amplifier as a phase detector, a PI controller as a loop filter, and the system is digital (on an FPGA embedded system).

A step response test is an easy way to determine the bandwidth. Sum a small step  into the control voltage of your oscillator (VCO or NCO), and measure the 90% to 10% fall time of the corrected response at the output of the loop filter as shown in this block diagram. Note that the loop will respond in such a way to completely cancel the injected offset, but can only respond at a rate within its bandwidth. Given this is in an FPGA implementation, this would be an easy test feature to have permanently included, especially valuable if this is a mixed signal system with external analog components that may vary from unit to unit or over temperature.

The bandwidth is given by the following approximation (this is accurate for a first order loop but will provide a reasonable estimate for higher order loops as well):
$$BW = \frac{0.35}{t} \space\space\space\space \text{(first order system)}$$
Where $BW$ is the 3 dB bandwidth in Hz and $t$ is the 10%/90% rise or fall time. This equation comes from solving the first order response equation given by $y(t) = e^{-t/\tau}$ where $\tau$ is the reciprocal of the bandwidth in radians/sec. For second order loops with a typical damping factor of 0.7 this relationship is closer to:
$$BW = \frac{0.33}{t}\space\space\space\space \text{(second order system, damping factor = 0.7)}$$
You could use this test for further loop forensics such as damping factor, natural frequency and a more accurate estimate of the bandwidth through observation of the overshoot and ringing in the corrected response.
The step should be as small as possible to stay within the tracking range of the loop, but large enough to make a reasonable measurement of the response time.
Below demonstrates the accuracy of using this first order approximation on a second order system, showing the step response and frequency response of a generalized 2nd order system with a natural frequency of 1KHz where the damping factor was varied between .5 and 1.0. The dashed red lines indicate the 10% to 90% crossings for rise time and the -3 dB BW crossing for the frequency response.
The following table shows the resulting accuracy of the approximation versus actual 3 dB bandwidth for each of these cases.


Can't get consistent magnitude response in MATLAB I have designed a notch filter with MATLAB filter designer app. Then I generated the filter design function which gives me the filter object:

function Hd = filter_notch

Fs=500; Fnotch=50; BW=1; Apass=1;

[b, a] = iirnotch(Fnotch/(Fs/2), BW/(Fs/2), Apass);
Hd     = dfilt.df2(b, a);

Now, I use freqz() to plot the response. However I get 2 different responses.

Mehtod 1. Direct plot with freqz(filter_object)

>> freqz(filter_notch);

Method 2. Manual plot with h=freqz(filter_object, Physical_frequencies, F_sample)

>> fs=500;
>> f_range=linspace(0,fs/2,2048);
>> h=freqz(filter_notch,f_range,fs);
>> hdb=20*log10(abs(h));
>> plot(f_range,hdb)

you can see the notch attenuations are different with above 2 methods, roughly -25dB vs -15dB. Why is that?

From MATLAB documents:

 h = freqz(___,f,fs) returns the frequency response vector, h, at the
 physical frequencies supplied in f.

As already mentioned by @Fat32 and @Florent Ecochard what your observing is a resolution issue with the produced plots. By default (e.g. using freqz without additional input arguments) 8192 points are used to plot the filter response; using less points, as done in you method 2, results in the observed behavior.

See the example code and figures below for an illustration:

%% Define Notch Filter
%
Fs     = 500;
Fnotch =  50;
BW     =   1;
Apass  =   1;

[b, a] = iirnotch(Fnotch/(Fs/2), BW/(Fs/2), Apass);
Hd     = dfilt.df2(b, a);

%% Plot Response

% method 1
% -> FREQZ default settings
%
[H1,F1] = freqz(Hd);
F1      = F1 * (Fs/(2*pi));  % convert to Hz / account for frequency normalization

% method 2
% -> Use frequency axis of 2048 points.
%
fs      = 500;
f_range = linspace(0,fs/2,2048);

[H2,F2] = freqz(Hd,f_range,fs);

% method 3
% -> Use frequency axis of 8192 points (same number as used in example 1).
%
fs      = 500;
f_range = linspace(0,fs/2,8192);

[H3,F3] = freqz(Hd,f_range,fs);

figure;
plot(F1,20*log10(abs(H1)),...
     F2,20*log10(abs(H2)),...
     F3,20*log10(abs(H3)) ...
    )
xlabel('f / Hz')
ylabel('|Hd| / dB')
legend('method 1 (freqz default)','method 2 (2048 points)','method 3 (8192 points)')

This then results in the following plots:


Image Projection with SIFT Recently, I am thinking about a problem regarding image projection. I have a larger picture as bellow.

I want to find the transformation matrix, x'=Hx. When the above image pass the transformation matrix H, I will get another smaller picture as bellow.

I have tried to use SIFT to find four major points and use these points to implement DLT algorithm. However, after executing SIFT, I find that there are many matching error as bellow. The distRatio is 0.5

I have thought for almost two weeks. Is there any different ideas which can be use to solve this problem. Thank you very much.

The logo you are looking at is symmetric. For example, you could rotate it by 120 degrees and get the same pattern. If you take a look, many of the "errors" are not really errors. The same pattern appears in your image many times. It is called self-similarity.

You could do several things:

Change descriptors to avoid confusion between features that are similar in an image.
Instead of finding one match for keypoint, find 3 matches and test multiple combinations by fitting a transformation and checking which fits better.
Limit the transformation by rotation angle, mirror, etc. For example, angle must not be larger than 60 degrees


Delays and alignment issues: how to compute the an equalizer delay? For example, if an equalizer model experiences a delay due to its processing, the equalizer may misinterpret the boundaries between the equalized symbols it outputs and, therefore, return meaningless results.

How to align the output data by finding the right delay?

Here a example of a simple equalizer implemented in matlab:

clc; clear all; N = 1e3; MQAM = 2;

d = randi([0 1],N,1);
x = qammod(d,MQAM,'gray','InputType','bit');
h = [0.1 0.5, 0.1 0.2];
r = awgn(complex(filter(h, 1, x), zeros(size(x))), 10, 'measured');

eq = comm.DecisionFeedbackEqualizer;
eq.Constellation = qammod(0:1, 2);
eq.StepSize = maxstep(eq,r);
eq.ReferenceTap = 2;

[y,e,w] = eq(r, x);

dec = qamdemod(y,MQAM,'gray','OutputType','bit');

errors = sum( dec~= d)

How to have dec and d align, in order word what is the delay ?

There are several ways to handle this, and specific details would depend on the equalization approach- but in post-processing applications with least squared equalizer solutions I either do a cross-correlation of the received and transmitted signal to understand both the overall delay as well as the delay spread, since using an equalizer span that exceeds the delay spread results in noise enhancement. Another approach I take is to just use a much longer span in the equalizer and then from the equalizer solution I can see where the dominant taps are, and from that understand both the delay and delay spread and narrow the equalizer span accordingly.

See example code implementations here where the equalizer span was set to assume dominant taps in the center (so was zero padded before and after), and this post demonstrates further the effect of delay:

Compensating Loudspeaker frequency response in an audio signal

And here where 0 delay was assumed since the equalizer approach was used to determine the delay specifically:

How determine the delay in my signal practically


channel impulse and frequency response: H? The received signal of MIMO system is given as:
$$ y= H \cdot x +n, (*)$$
where $H$ is a channel matrix.

$h(\tau_l)$ is a channel impulse response. It is determined as a matrix $\in \mathbb{C}^{Nr \times Nt}$, $\tau_l$ is delay.

My question: $H$ in (*) is a frequency response, fourier transform of   $h(\tau_l)$, isnt't?

No, $H_{rx,tx}$ is the channel matrix with $(r,t)^{th}$ element representing channels attenuation and phase rotation on transmit symbols $\vec{x}_{tx}$, from $t^{th}$ transmit antenna to $r^{th}$ receive antenna. The equation is neither in time-domain nor in frequency domain. It is on transmit symbols. For example, in OFDM systems, you get $y_{rx}$ (received symbols) only after FFT operation at the receiver, so you can understand the equation as in frequency domain. But for a Non-OFDM based system, like CDMA, you get received symbols $y_{rx}$ after demodulation with matched-filter and that would be in time-domain.

You can understand this equation as a relationship between $rx$ number of received symbols $y_{rx}$ and $tx$ number of transmitted symbols $x_{tx}$ in MIMO scenario where $n_{rx}$ is the AWGN for each of the $rx$ receive chain and $H_{rx,tx}$ is the channel matrix affecting the symbol transmitted from $t^{th}$ antenna and received at $r^{th}$ antenna.
Since symbols are represented as vector points complex plane, hence each channel can only affect each symbol by changing its magnitude and giving a random phase rotation. It doesnot matter time-domain or frequency-domain. That depends on your Receive chain implementation of getting $\vec{y_{rx}}$.


Getting error in fit I was wondering if anyone had run into the problem of trying to estimate errors in their signal processing on spectroscopic data. I know that many people use spectroscopic techniques to estimate concentrations of different materials, but I would like to know how accurate these measurements can be  (i.e. I would like to get an output something like "species A is 20% +-1% of the sample given").  Further, I would like to know how to deal with very extreme cases, where several different types of materials are present, and may have observables which fall directly on top of each other.

A simple example may be the following:

You can see there are two species being fit to the data. If the areas are then transformed by calculation into percentages, (i.e. the sample taken is 48% B, 52% A) how can we be sure of this and how accurate is the fit?   I know this will be dependent on the accuracy of the estimate of the positions of these peaks that are (perhaps) given by the user, so I am interested in a method that takes a known error in the parameters (say +-15 on the x axis for error in the peak center position, +-10 error in the width of the peak, etc.).

I suspect that the errors will become large when the observables overlap (i.e. the peak centers for two fitting functions are the same).

In addition, it is possible that these spectra have a large background, which may also have error, affecting all of the other species and their errors.  I am not certain if this background would be treated differently than all of the other species, or if it could be treated within the same algorithm as all of the other species.

To illustrate my point further, here is an image of a spectroscopic measurement of several different materials:

On the top in red is the raw data given by optical absorption (the measured data), while the black is a calculated background, and the blue is a calculated sum of all of the species including the background.

On the bottom, the calculated background is subtracted (blue and red lines now DO NOT include the calculated background), while the several different colored lines below are each individual species being summed up to create the blue line.

These are the calculated measurements I am interested in estimating the error in.

As you can see, the error is enormously large in this example for most of the calculated measurements. Each species may or may not have several 'peaks' associated with it, which can be illustrated by the bolded yellow calculated line.  In addition, you can see that several of the calculated peak centers fall around the same place, so this will likely reduce the certainty that the measurements are correct even if the calculated line falls directly upon the raw data.

I have calculated the mean squared displacement as a quick estimate of how good the fit is, but I know that this doesn't do anything to address any larger concerns of the actual calculated measurement uncertainty. The most I have really done in statistics is standard deviation and calculations dealing with several measurements, but this is quite different, since it deals with how sure you can be with only one measurement, not seeing differences in multiple measurements. Is there a way of using confidence intervals and confidence levels on non independent and identically distributed random variables?  (Again, I am very new to statistics and have never taken a course on it, so I apologize if this is elementary or trivial)

Well, without knowing the specifics of your algorithm, I assume that you basically have multiple functions (e.g. Gaussian or Lorentzian peaks, each with position and FWHM parameters, additionally maybe some polynomials for backgrounds etc.) and you sum them all up into one big "fit function" that you hand over to an optimization algorithm that jiggles the parameters around until you found a good fit in the sense of non-linear least squares.

What you can do then is computing the asymptotic standard errors. For this, you have to (most likely numerically) compute the Jacobi Matrix of your fit function (matrix of partial derivatives for each parameter (columns) at each measurement point (rows)). From this you can compute the variance-covariance matrix and then compute the errors for each parameter. This is an approach very often found in curve fitting software. In fact, there is a quite comprehensive explanation here at OriginLabs (see the section on Parameter Standard Errors). Also a more detailed overview can be found here at arxiv.

Other methods are discussed in this paper. They certainly are appropriate and maybe even mandatory in some cases, but I find them a bit too complicated for most "real world" applications.

If I were you, I would go with the asymptotic standard errors for a first start and see, if they fit the needs. It is an easy to implement strategy that is, as already said, also widely found in software packages.


Exponential moving average cut-off frequency I am trying to implement low pass filter from this example.

What is the cut-off frequency for this type of filter? Is it $F_s \frac{1-\alpha}{2\pi\alpha}$, where $F_s$ is sampling frequency?

The discrete recurrence relation given on the linked page is
$y[n] = (1-\alpha)y[n-1] + \alpha x[n]$
with $x[n]$ being the input samples and $y[n]$ being the output samples.
So taking the Z-transform
$Y(z) = (1-\alpha)z^{-1}Y(z) + \alpha X(z)$
$H(z) = \dfrac{Y(z)}{X(z)}= \dfrac{\alpha}{1-(1-\alpha)z^{-1}}$
To find the -3 dB corner of the filter, in terms of the normalized digital radian frequency, $\Omega$ in the range $[0, \pi]$, one solves
$\left|H(z=e^{j\Omega_{3dB}})\right|^2 = \dfrac{1}{2} = \left|\dfrac{\alpha}{1-(1-\alpha)e^{-j\Omega_{3dB}}}\right|^2$
$ = \dfrac{\alpha^2}{\left|1-(1-\alpha)\cos(-\Omega_{3dB}) -j (1-\alpha)\sin(-\Omega_{3dB})\right|^2}$
which results in this equation
$\left[1-(1-\alpha)\cos(\Omega_{3dB})\right]^2+\left[(1-\alpha)\sin(\Omega_{3dB})\right]^2 = 2\alpha^2$
After some more algebra and trigonometry, I come up with
$\Omega_{3dB} = \cos^{-1}\left[1-\dfrac{\alpha^2}{2(1-\alpha)}\right]$
and since
$f_{3dB} = \dfrac{\Omega_{3dB}}{\pi} \cdot \dfrac{F_s}{2}$
I come up with
$f_{3dB} = \dfrac{F_s}{2\pi}\cos^{-1}\left[1-\dfrac{\alpha^2}{2(1-\alpha)}\right]$
Here's some Octave/MatLab code so you can verify the result on a plot; just change your Fs and alpha as required:
Fs = 10000.0
alpha = 0.01

b = [alpha];
a = [1 -(1-alpha)];

[H, W] = freqz(b, a, 1024, Fs);

f3db = Fs/(2*pi)*acos(1-(alpha^2)/(2*(1-alpha)))

plot(W, 20*log10(abs(H)), [f3db, f3db], [-40, 0]);
grid on

Update in response to comment:
The graph is a plot of the filter's magnitude response in dB (on the y-axis) vs. input frequency in Hz (on the x-axis).  The vertical line designates the -3 dB corner frequency.
Since you want to specify your -3 dB corner frequency to find the value for $\alpha$, let's start with this equation from above:
$\left[1-(1-\alpha)\cos(\Omega_{3dB})\right]^2+\left[(1-\alpha)\sin(\Omega_{3dB})\right]^2 = 2\alpha^2$
and after some algebra and trigonometry, one can get a quadratic equation in $\alpha$
$\alpha^2 +2(1- \cos(\Omega_{3dB}))\alpha - 2(1- \cos(\Omega_{3dB})) = 0$
which has solutions
$\alpha = \cos(\Omega_{3db}) - 1 \pm \sqrt{\cos^2(\Omega_{3dB}) -4\cos(\Omega_{3dB}) +3}$
of which only the $+$ answer from the $\pm$ can yield a positive answer for $\alpha$.
Using the above solution for $\alpha$ and the relation
$\Omega_{3dB} = \dfrac{\pi}{F_s/2}f_{3dB}$
one can find $\alpha$ given $f_{3dB}$ and $F_s$.
Here is an updated Octave script:
Fs = 40000
f3db = 1

format long;
omega3db = f3db * pi/(Fs/2)

alpha = cos(omega3db) - 1 + sqrt(cos(omega3db).^2 - 4*cos(omega3db) + 3)

b = [alpha];
a = [1 -(1-alpha)];

[H, W] = freqz(b, a, 32768, Fs);

figure(1);
plot(W, 20*log10(abs(H)), [f3db, f3db], [-40, 0]);
xlabel('Frequency (Hz)');
ylabel('Magnitude Response (dB)');
title('EWMA Filter Frequency Response');
grid on;

W2 = [0:75] * pi/(Fs/2); % 0 to 75 Hz
H2 = freqz(b, a, W2);
W2 = W2 / (pi/(Fs/2));

figure(2);
plot(W2, 20*log10(abs(H2)), [f3db, f3db], [-20, 0]);
xlabel('Frequency (Hz)');
ylabel('Magnitude Response (dB)');
title('EWMA Filter Frequency Response near DC');
grid on;

Additional information in response to Sep 11 '20 comment
For a stable, low pass filter, the pole, $(1-\alpha)$, must be inside the unit circle of the z-plane and be greater than $0$ on the real axis.  So we must satisfy
$$0 < 1 - \alpha < 1$$
thus
$$0 < \alpha < 1$$
These are bounds on $\alpha$ (and not a "definition") that come from requirements for stability and a filter passband centered at DC.
However, not all valid values of $\alpha$ will yield a -3 dB frequency that exists in the interval from DC to the Nyquist frequency, $\Omega \in [0, \pi]$ or equivalently $f \in \left[0, \dfrac{F_s}{2}\right]$.
When the pole, $(1-\alpha)$, is very close to the origin of the z-plane, the filter response will be so gradual or flat, that it will never drop to -3 dB or below.  Specifically, when
$$\Omega_{3dB} = \cos^{-1}\left[1-\dfrac{\alpha^2}{2(1-\alpha)}\right]$$
is undefined, due to
$$\left|1-\dfrac{\alpha^2}{2(1-\alpha)}\right| > 1$$
the filter response is so gradual, it will not have a drop down to -3 dB or below before the Nyquist frequency.  The $\alpha$ at which this happens is given by the positive solution to
$$ \alpha^2+4\alpha-4 =0$$
which is
$$\alpha = -2 + 2\sqrt{2} \approx 0.82843$$
EWMA filters with an $\alpha$ that is this high, have a pole so close to the origin of the z-plane, that the filter response is so gradual, there is no frequency attenuated by 3 dB or more.
The Wikipedia derivation of the 3 dB cutoff point is flawed and comes to an incorrect/inexact result.  Specifically it directly applies a continuous time, Laplace s-domain relationship:
$$f_c = \dfrac{1}{2\pi RC}  \quad \text{(Don't use this!)}$$
into a derivation of a result for a discrete time, z-domain, filter.  That step invalidates the derivation of an exact result in the Wikipedia article.
As a consequence, the result of the derivation in the Wikipedia article is only an approximation and (subjectively) only good for cutoff frequencies up to about $0.1 F_s$.


Speech Recognition - Subtitles I am looking for ANY advice anyone can give me regarding on this problem.

I'm looking to start a project on subtitles, basically, if we take a Youtube video of a breaking news story, subtitles are generated automatically without the need of physically inputting what the speaker is saying. Now, I understand the complexity of this problem, and, I am aware of the amount of training samples I would need in order to carry out such a project.

The question is this: I have had experience in using DTW and HMM in order to train and recognise simple words "Yes", "No" and gained a lot of experience from this and as part of my previous findings I found that the DTW algorithm is slow and therefore would not be efficient for this. For this project, I imagine real-time processing is involved. Would it therefore be a wise thing to implement a Hidden Markov Model for such a problem, or, would another model/algorithm be more suitable?

Here is my basic structure:

I will upload a 3 minute video to youtube, I will then extract the raw audio from this video and then extract the features (MFCC) which will then be used in order to train the model (HMM) which will /hopefully/ identify the spoken word.

Could anyone recommend any papers that could potentially help me with tackling such a problem?

Thank you for taking time to read this. If anything is not understood, please let me know.

The amount of knowledge necessary to develop such a large scale multi-language, speaker-independent, large-vocabulary speech recognition system is spread well over hundreds of papers; and each individual brick of the system (say the feature extraction front-end, the FST decoding library, the language model store) is developed by world-class expert in this specific field. Some of Google's research scientists have published papers on these topics, this might be a good place to start if you want to understand what's behind the scenes... But maybe  you should start with a smaller goal - like building a system that has a vocabulary of a dozen different words from the same speaker - using a library like HTK or Sphinx (using the "canonical" approach with the hierarchy GMM->phone->word->language model), and once you have all the basics right, explore deeper in one direction (try more complex feature extraction methods to handle noise, grow the vocabulary and the complexity of the grammar, try speaker adaptation).


MATLAB flanger...zipper noise, crackling, help with interpolation values? i'm writing a batch script in MATLAB to flange.

i'm getting some periodic noise, and i'm having trouble tracking down the source. it's definitely in the algorithm, because soundsc() will scale values larger than 1. and theres no memory or computing time issue, because it's running as a batch process.

any mp3 or wave file should do the trick

current script, with interpolation (sounds worse so far than non-interpolating, because i am probably using the wrong X for the interpolation.)

function result = dLinTerp(x1,x2,y1,y2,x)

% calculate decimal position of x
dx = (x - x1)/(x2 - x1);

% use weighted sum method of interpolating
result = dx*y2 + (1-dx)*y1;

and

    cur_frac=cur_sin*max_samp_delay;

    if (readIdx + cur_frac >= (max_samp_delay))
        fetchIdx = cur_frac - (max_samp_delay - readIdx) + 1;
    else
        fetchIdx = readIdx + cur_frac;
    end

    if (floor(fetchIdx) = max_samp_delay)
        fetchNext = 1;
    else
        fetchNext = ceil(fetchIdx);
    end
     integ = fix(cur_frac);
     frac = abs(cur_frac - integ);
    interp = dLinTerp(0.0, 1.0, delayBuff(fetchPrev), delayBuff(fetchNext),frac);
    out(i) = (amp*input(i)) + (amp*interp);

here is the original script, no interpolation

%flanger

[input,fs] = audioread('pianoriff_s.mp3');
input = input(:,1);
len = length(input);
out = zeros(len,1);

max_time_delay=0.007; % 7ms max delay in seconds
rate=0.2; %rate of flange in Hz
max_samp_delay= floor(max_time_delay*fs);
delayBuff = zeros(max_samp_delay,1);
amp=0.7;
readIdx = 1;

for i = 1:len,
    delayBuff(readIdx) = input(i);

    sin_ref = sin(2*pi*i*(rate/fs));
    cur_sin=abs(sin_ref); %abs of current sin val 0-1
    cur_delay=ceil(cur_sin*max_samp_delay);

    if (readIdx + cur_delay >= max_samp_delay)
        fetchIdx = cur_delay - (max_samp_delay - readIdx) + 1;
    else
        fetchIdx = readIdx + cur_delay;
    end

    out(i) = (amp*input(i)) + (amp*(delayBuff(fetchIdx)));

    if (readIdx+1 >= max_samp_delay)
        readIdx = 1;
    else
        readIdx = readIdx + 1;
    end
end

soundsc(out,fs);

ok so it turns out setting indices is tricky in MATLAB (who knew :p)

first, a less needlessly complicated interpolator:

function result = dLinTerp(y1,y2,x)

% use weighted sum method of interpolating
result = x*y2 + (1-x)*y1;

and the flanger

% flanger
clear all;
% Creates a single FIR delay with the delay time oscillating from
% Either 0-3 ms or 0-15 ms at 0.1 - 5 Hz
[input,fs] = audioread('pianoriff_s.mp3');
input = input(:,1);
len = length(input);
out = zeros(len,1);

% parameters to vary the effect %
max_time_delay=0.005; % 7ms max delay in seconds
rate=0.3; %rate of flange in Hz

% convert delay in ms to max delay in samples
max_samp_delay= floor(max_time_delay*fs);
delayBuff = zeros(max_samp_delay,1);
amp=0.7;
readIdx = 1;

fetchPrev = 1;

% delaybase = .010;
% new_d_base = delaybase*fs;
% d_base = LIN_INTERP(0,1,0, new_d_base,frac);
%
% dp = (float)(delay_pos - d_base) - (delay_depth * law);
% %     // Get the integer part
% dp_idx = f_round(dp - 0.5f);
% %     // Get the fractional part
% dp_frac = dp - dp_idx;

% for each sample
for i = 1:len,
    delayBuff(readIdx) = input(i);
    sin_ref = 0.97*sin(2*pi*i*(rate/fs));
    sin_next = 0.97*sin(2*pi*(i+1)*(rate/fs));
    cur_sin=abs(sin_ref); %abs of current sin val 0-1
    cur_delay=ceil(cur_sin*max_samp_delay);

    cur_frac=cur_sin*max_samp_delay;

    if (readIdx + cur_frac >= (max_samp_delay))
        fetchIdx = cur_frac - (max_samp_delay - readIdx) + 1;
    else
        fetchIdx = readIdx + cur_frac;
    end

    if (floor(fetchIdx) = max_samp_delay)
        fetchNext = 1;
    else
        fetchNext = ceil(fetchIdx);
    end

    integ = fix(cur_frac);
    frac = abs(cur_frac - integ);
    interp = dLinTerp(delayBuff(fetchPrev), delayBuff(fetchNext),frac);
    out(i) = (amp*input(i)) + (amp*interp);

    if (readIdx+1 >= max_samp_delay)
        readIdx = 1;
    else
        readIdx = readIdx + 1;
    end

end
% write output

soundsc(out,fs);


Low-pass vs median filter delay I need to design a filter to reduce noise from accelerometer readings and use this data to a control system. One of the problems of a simple low-pass filter (like a single pole IIR filter) is the delay introduced by the filter. As an alternative I saw many people using a kalman filter or a median filter. I'm very interested in this median filter due to its simple formulation. But then I realize that it will always have some samples lag in the output. For the filter being able to output the present term, it should be able to read N/2 (or (N-1)/2) terms from the future which is obviously not possible in real time systems. So, at best, its output would be N/2 samples late.

The only difference I see in terms of delay between these two filters is that the low-pass may have different delays for different frequencies input while the median-filter has a constant sample lag. The larger the window size is(which is also related to how smooth our output will be), larger will be this lag. So it might be interesting to reconsider a single pole filter.

Is that reasoning correct?

i would say that your reasoning is a teeny bit flawed.

the delay of an LTI filter is a deterministic function (of frequency) that is unambiguous.  a sinusoid of some known frequency goes in and the sinusoid that comes out will be delayed by a predetermined amount of time.

a median filter is not linear.  and, if a sinusoid goes in, it will not be a sinusoid coming out.  let's say that your causal median filter is outputting the median value of the current input sample and the most recent past $N-1$ samples (the "formulation" may seem "simple", but you're gonna need to repeatedly sort $N$ samples to pick out your median value), you can come up with an expected value of the delay and, since the sample that is eventually chosen to be the median value can be any one those $N$ samples with equal likelihood, then the expected value of the delay is $\frac{N-1}{2}$.


How do I implement cross-correlation to prove two audio files are similar? I have to do cross correlation of two audio file to prove they are similar.
I have taken the FFT of the two audio files and have their power spectrum values in separate arrays.

How should I proceed further to cross-correlate them and prove that they're similar? Is there a better way to do it? Any basic ideas will be helpful for me to learn and apply it.

Cross-correlation and convolution are closely related.  In short, to do convolution with FFTs, you

zero-pad the input signals (add zeros to the end so that at least half of the wave is "blank")
take the FFT of both signals
multiply the results together (element-wise multiplication)
do the inverse FFT

conv(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros))

You need to do the zero-padding because the FFT method is actually circular cross-correlation, meaning the signal wraps around at the ends.  So you add enough zeros to get rid of the overlap, to simulate a signal that is zero out to infinity.

To get cross-correlation instead of convolution, you either need to time-reverse one of the signals before doing the FFT, or take the complex conjugate of one of the signals after the FFT:

corr(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros[reversed]))
corr(a, b) = ifft(fft(a_and_zeros) * conj(fft(b_and_zeros)))

whichever is easier with your hardware/software.  For autocorrelation (cross-correlation of a signal with itself), it's better to do the complex conjugate, because then you only need to calculate the FFT once.

If the signals are real, you can use real FFTs (RFFT/IRFFT) and save half your computation time by only calculating half of the spectrum.

Also you can save computation time by padding to a larger size that the FFT is optimized for (such as a 5-smooth number for FFTPACK, a ~13-smooth number for FFTW, or a power of 2 for a simple hardware implementation).

Here's an example in Python of FFT correlation compared with brute-force correlation:  https://stackoverflow.com/a/1768140/125507

This will give you the cross-correlation function, which is a measure of similarity vs offset.  To get the offset at which the waves are "lined up" with each other, there will be a peak in the correlation function:

The x value of the peak is the offset, which could be negative or positive.

I've only seen this used to find the offset between two waves.  You can get a more precise estimate of the offset (better than the resolution of your samples) by using parabolic/quadratic interpolation on the peak.

To get a similarity value between -1 and 1 (a negative value indicating one of the signals decreases as the other increases) you'd need to scale the amplitude according to the length of the inputs, length of the FFT, your particular FFT implementation's scaling, etc.  The autocorrelation of a wave with itself will give you the value of the maximum possible match.

Note that this will only work on waves that have the same shape.  If they've been sampled on different hardware or have some noise added, but otherwise still have the same shape, this comparison will work, but if the wave shape has been changed by filtering or phase shifts, they may sound the same, but won't correlate as well.


Why does the Fourier Transform of the impulse look so different from the Fourier Transform of the impulse train? The fourier transform of the impulse functions is:

$$ \delta(t) \longleftrightarrow 1$$

The shifted delta:

$$ \delta(t-nT) \longleftrightarrow e^{-j \Omega nT}$$

But the fourier transform of the impulse train is:

$$\sum_{n=-\infty}^{\infty} \delta(t-nT) \longleftrightarrow \frac{2\pi}{T} \sum_{k=-\infty}^{\infty} \delta \left(\Omega - k\tfrac{2\pi}{T} \right)$$

How could I go about proving the impulse train fourier transform using the definition of the fourier transform of the shifted delta? It seems confusing... Thank you!

This answer shows that $\sum\limits_{n=-\infty}^{\infty}
  e^{j  n \Omega T}$ is actually the Fourier Series of $\frac{2\pi}{T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\tfrac{2\pi}{T})$.

\begin{equation}
 \delta(t-nT) \frac{2\pi}{T} e^{j  n \Omega T}
\end{equation}
which means
\begin{equation}
 \sum\limits_{n = -\infty}^{\infty} \delta(t-nT) \longleftrightarrow \sum\limits_{n = -\infty}^{\infty} e^{j  n \Omega T}
\end{equation}
Now consider the exponential Fourier series of $\frac{2\pi}{ T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)$, because this function is periodic, so we can say
\begin{equation}
 \frac{2\pi}{ T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)
 =
 \sum\limits_{n=-\infty}^{\infty}
 c_n e^{j  n \Omega T}
\end{equation}
where
\begin{equation}
 c_n = \frac{1}{\Omega_s}\int\limits_{-\frac{\Omega_s}{2}}^{\frac{\Omega_s}{2}}
 \frac{2\pi}{ T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)
 e^{j  n \Omega T}
 d (\Omega T)
\end{equation}
and $\Omega_s = \frac{2\pi}{T}$ by definition.

Take $\alpha = \Omega T$,
\begin{equation}
 c_n = \frac{1}{\Omega_s}
 \frac{2\pi}{ T} \sum_{k=-\infty}^{\infty}
 \int\limits_{-\frac{T\Omega_s}{2}}^{\frac{T\Omega_s}{2}}\delta(\frac{\alpha
}{T} - k\Omega_s)
 e^{-j  n \alpha}
 d \alpha
\end{equation}
The term $\frac{1}{\Omega_s}
 \frac{2\pi}{T}  = 1$ and take now
\begin{equation}
 \beta = \frac{\alpha}{T}
\end{equation}
\begin{equation}
 c_n =
  \sum_{k=-\infty}^{\infty}
 \int\limits_{-\frac{\Omega_s}{2}}^{\frac{\Omega_s}{2}}\delta(\beta - k\Omega_s)
 e^{-j  n \beta T}
 d \beta
 \tag{1}
\end{equation}
We know that
\begin{equation}
 \int_{a}^b \delta(x - x_0) f(x) \ dx
 =
 f(x_0)
 \quad
 \text{if } a < x_0 < b
\end{equation}
and zero elsewhere. So the integral in $(1)$ contains one value within $[-\frac{\Omega_s}{2},\frac{\Omega_s}{2}]$ i.e. for $k = 0$, so
\begin{equation}
 c_n =
  \sum_{k=-\infty}^{\infty}
 \int\limits_{-\frac{\Omega_s}{2}}^{\frac{\Omega_s}{2}}\delta(\beta - k\Omega_s)
 e^{-j  n \beta T}
 d \beta
=
e^{-j  n (0) T}
=
1
\end{equation}
So $c_n = 1$ for all $n$. We have just proved that
\begin{equation}
 \frac{2\pi}{ T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)
 =
 \sum\limits_{n=-\infty}^{\infty}
  e^{j  n \Omega T}
\end{equation}
So whether you say $\sum\limits_{n=-\infty}^{\infty}
  e^{j  n \Omega T}$ or $\frac{2\pi}{ T} \sum_{k=-\infty}^{\infty} \delta(\Omega - k\Omega_s)$, it is the same thing.


Confusion in linearity? I am studying linear systems from the book "Signals and Systems Laboratory with MATLAB". I have performed an example on MATLAB myself in this regard.
Code is as follow:

clc
clear all
close all
t=-3:.1:3;
x1=sin(t);
x2=cos(t);
a1=2;
a2=3;
z=a1*x1+a2*x2;
y1=2*z+2*cos(pi*t/3)
plot(t,y1);
ylim([-1 11]);
figure;
z1=2*x1+2*cos(pi*t/3);
z2=2*x2+2*cos(pi*t/3)
y2=a1*z1+a2*z2;
plot(t,y2);
ylim([-1 19]);

When I run this program I get two  figures containing two same shapes that are curve-like but amplitude/size/height is not same.

I am confused if my this system defined by equation y=2*z*cos(pi*t/3) used in program above is linear or not? If still yes, then why (as we have size difference in both curves)?

To expand on Justme's reply:

From the looks of your code you seem to have confused what constitutes an LTI system.

In your code the LTI system that I see (y1=2*z) scales an undelayed input signal by 2. The additional 2*cos(pi/t) is NOT part of the system, it is in effect just another sinusoid that you have added to your signal - thus giving it an additional frequency component.

An LTI system can only ever do two things to an input complex signal - it can scale the magnitude and it can modify the phase of frequency components of the signal. It cannot modify the fundamental frequencies of the sinusoidal components of a signal. (NB: It can however create resonances that self-oscillate - this is however cause by recursion in the system)

Returning to your code once again, we see that you have checked the linearity property, but have also included a scaling of the additional sinusoid that you have misattributed as part of the system. If you remove this sinusoid in your system function, you will find that the linearity property holds for your system

Edit: Fixed some technical wording issues


Understanding Pitch Detection with Autocorrelation Methods I've been reading through A Smarter Way to Find Pitch, describing its pitch detection algorithm using autocorrelation.

I've been having trouble understanding the accuracy claims. It says:

  MPM runs in real time with a standard 44.1 kHz sampling rate. It operates without using low-pass filtering so
  it can work on sound with high harmonic frequencies such
  as a violin and it can display pitch changes of one cent reliably.

I cannot fathom how detecting changes of one cent is possible. Considering a note played at 441Hz, sampled at 44.1kHz, one would measure a first peak at tao=100 samples (delay space).

The next closest note which could be measured, by my understanding, would have a peak at the 101st sample. This would be a frequency of 437Hz, or about 16 cents lower than the first note. What am I missing here?

For higher resolution periodicity estimation, given properly band-limited data, the data can be upsampled (interpolated) before autocorrelation, and/or, the autocorrelation results can be interpolated (Sinc kernel or some order of polynomial).

This, of course, assumes the signal has a high enough S/N, and that the periodicity itself is unmodulated in timbre(harmonic content), phase, frequency, or amplitude, which isn’t always true of real-world pitched sounds.  Thus a larger number of periods may be needed for a some desired statistical error bound on the pitch estimation.


How to programmatically find a sharp change in slope First, apologies if I have the wrong group or this question is far too easy for this group.  I am (as you'll see) a newbie. Please point me to another, more appropriate group OR tell me how to fix the question to be more appropriate for this group.

I have some data that represents a thermometer that is put in a fridge, and then a freezer, and then finally, back out of the freezer to the fridge.  The data might look like this:

The horizontal axis is time in seconds and the vertical axis is Celsius. As you can see, it's quite obvious where the thermometer goes into the freezer and where it comes out.  One can with the naked eye come up with a pretty good  estimate as to where this happens.  I'd like to be able to have a computer program (C#) figure this out for me for other data sets that all have this general form.  One way  might be to look for when the data first reaches -10 and argue that the transition happened at about this time.  Another might be to look for a change in the 2nd derivative.  That's about the extent of my knowledge.  I can imagine there is a whole literature devoted to this problem and probably libraries of code that already do this. I'm looking for

1) Starter ideas on what to read to educate myself about this problem.

2) Example code in C# (or JAVA or C++) that does the automatic detection.

I would like to be able to feed a time series into an algorithm, and have an estimate of when the thermometer was put in the freezer (or taken out) WITH some sort of confidence interval. Example: "it went into the freezer at 18420 seconds +/- 100 seconds with 95% confidence".  I am willing to assume (for now) that the thermometer starts changing temperature as soon as it is moved from the fridge to the freezer (I'm still researching this!).

Those are the questions.  For those, who are curious for background information, please continue to read.  I'll add that the motivation is that we are trying to see if some thermometers in the fridge and freezer are in spec.  We have several thermometers in the fridge and freezer that we would like to know are in spec.  We have a highly accurate "reference" thermometer that we can move around with a robot to the locations where the fridge/freezer thermometers are.  As we move the reference thermometer to the various locations with the robot, we have the relevant fridge or freezer thermometer take a temperature and record it.  This we can do with very accurate timing.  The problem is that the reference thermometer just has a counter every 10 seconds.  Presumably the measurements are really 10 seconds apart but there is not absolute time measurements.  It's just a counter from 0.  You can see the graph of the reference thermometer in my attached picture.  It occurs to me that since I move the reference thermometer to the first freezer station, I should be able to use that point on the graph as a "pin" to connect the two timelines.  I'd further point out that I leave the reference thermometer at each station for several minutes, so there is some slop available.  Finally, for those wise folks that say, "why not have the reference thermometer communicate with the rest of the system via Bluetooth or similar so you don't have to try and match things up after the fact?", my reply is "already suggested (several times) to management without success" :)

This looks like an interesting problem.

First, if your timing is evenly spaced, your idea of having two "pin" (aka calibration) points should give you an adequate linear relation between your two time scales.  If you have more than two points then you can use linear regression to find a best fit.

Second, I think you will want to use the raw values and values of a calculated first derivative to classify parts of your graph.  The actual value of a second derivative doesn't seem to me to be particularly helpful.  Because of the noise you have, using immediate first difference calculations may lead to bad results.

I have a suggestion for a different approach which may work quite well.  Check out the difference technique in my blog article Exponential Smoothing with a Wrinkle.  You need not necessarily use exponential smoothing, boxcar smoothing or some other would also work.  The idea is to smooth in a forward direction and also smooth in a backward direction.  At every point, the forward smoothing gives you the "expected value" from the past and the backward smoothing gives you the "expected value" from the future.  Substracting the former from the latter gives you a first derivative like signal.  The coding is rather trivial so I will not provide it.  I suggest you give it a shot on your data and see what the graphs look like.  The smoothing will solve your noise issue.

I'm not sure the best way to calculate your probability values, but I think with a few graphs and calculations under your belt you will be able to provide an estimate that will be valid and fairly consistent across all your cases.


Can we have a Digital Anti Aliasing filter? I am working on a board that has no antialisaing filter at the input of the ADC. I have  option to I implement my own filter using RC + Opamp circuit. But is it also possible to implement Anti Aliasing filter after sampling by ADC and processing in Digital domain: a digital Anti aliasing filter?

I agree with pichenettes's answer but I would like to add that it is pretty common practice to use a simple inexpensive low-order analog anti-aliasing filter, and do the rest of the anti-aliasing filtering in the digital domain. This implies of course that you do not process at maximum sampling rate but that you downsample after the digital anti-aliasing filter. Summarizing:

Of course you need an analog anti-aliasing filter.
The analog filter can be kept very simple if you can downsample your signal. In this case you can do more aliasing removal in the digital domain (before downsampling).


Why call it discrete time signals as difference equations? There are lots of representations of discrete time signals. You represent it in block diagram or as difference equations.

We call it block diagrams as there are elementary structures in the diagrams.

Why do we call the equation part as difference equation of a signal?
$y[n] = x[n] + x[n-1]$

There is nothing getting subtracted in the above example representation, Then why call it difference equation??

A difference equation is a recurrence relation where the current element $y[n]$ of a sequence is related to its past values $y[n-k]$, $k>0$ (if we interpret the index $n$ as time index, which is not necessary). The equation you gave in your question is a special case, because there are no past terms $y[n-k]$ involved, so you could call it a zeroth order difference equation, but you're right that from a common sense point of view it is no difference equation at all.

However, what counts here is the fact that no past (or future) elements $y[n-k]$ are involved, and not that "nothing is getting subtracted". E.g., the following equation is a standard example of a difference equation:

$$y[n]=ay[n-1]+bx[n]\tag{1}$$

The subtraction takes place in the index $n$ ($y[n-1]$). But obviously you can substitute $n+1$ for $n$ yielding

$$y[n+1]=ay[n]+bx[n+1]\tag{2}$$

which is completely equivalent to $(1)$ with no subtraction taking place, but it's still a difference equation, because the difference between the indices remains unchanged.


A good ecc algorithm for random errors? I know I will have a 10% error rate, but this errors are truly random. I was using solomon reed codes, but if a single bit in the code or parity word is bad the whole word becomes useless. If I have eight bytes of data and eight bytes of parity codes then 5 bit error can make my whole message bad. Is there a better algorithm for cases like this?

Edit:

I am trying to send packets of 128bit or 16bytes. 1/2 is data 1/2 is parity. I am having trouble figuring out how LDPC would be more beneficial.

If I have the same 16bytes with 8 bytes of parity and 8 bytes of data how ldpc be better? For example suppose with my example I have 1 bit errors in 5 different bytes. This would make SR unable to decode my data. Would it be the same LDPC or would it be able to recover that. From would I understand from LDPC it would still be unable to decode it.

First, it is not true that an error in the parity portion of a Reed-Solomon code word renders the word useless.  RS codes can correct errors in their parity and their data.

Having said that, 10% errors is very high and I would not use a RS code in that situation.  In fact, in a situation with that many errors you will occasionally have > 50% errors just due to random distribution (random distribution is often "clumpy") of the errors.  No code you use will correct all of the errors all of the time.  Your system should be designed with that fact in mind.

I would use a big turbo code or LDPC code.  Either one will get you near optimal performance.  The code needs to be "big" in the sense of using large blocks to avoid the clumpiness problem noted above.  You will still get more errors on occasion, but (assuming a gaussian distribution) your tails will be MUCH smaller if you make your blocks large.

If your errors are bursty you can use interleaving to spread them out among multiple blocks.  If they aren't bursty, interleaving won't help any.


Why in MuSiC (DoA) don't we search for minima of the inverse (pseudo)spectrum? I have started working on a Direction-of-Arrival estimation problem and I am testing various algorithms. I came across Multiple Signal Classification (MuSiC) method. In all textbooks I have looked at, the final step of the method is to search for maximums in the pseudo-spectrum given by
$$ P \left(\theta \right) = \frac{1}{a^{H}\left(\theta\right) V_{n} V_{n}^{H} a\left(\theta\right)}$$
where $P\left(\theta\right)$ is the spatial pseudo-spectrum, $a\left(\theta\right)$ is the array manifold of the array under consideration, $[ \cdot]^{H}$ denotes Hermitian transposition and $V_{n}$ is the matrix whose columns form the noise subspace (the noise eigenvectors).
Now, my question is, why should we look for maxima in $P\left(\theta\right)$ an not for minima in $\frac{1}{P\left(\theta\right)}$? Is there something I am missing here? I believe that the implementation of the algorithm could potentially benefit from avoiding this division.
Any insights and/or hints are most welcome. Thanks in advance.

Let's define the null-spectrum $Q(\theta)$ and pseudo-spectrum $P(\theta)$ as
$$Q(\theta) = a(\theta)^H{v_n}{v_n^H}a(\theta)$$
$$P(\theta) = \frac{1}{Q(\theta)}$$
The choice to find the maxima in $P(\theta)$ comes down to personal preference. It is more intuitive to associate a peak in a signal to some phenomenon being observed. This is especially so when reading a text, where the author is trying to teach you how the algorithm works.
You can very well search for the minima in $Q(\theta)$ instead of the maxima in $P(\theta)$ and get the same answer as can be seen in the example below

Nothing is stopping you from using the minima. In general there is no significant advantage to using either except for avoiding division operations and potential divide-by-zero, which usually fine in real systems where some kind of noise is present. Still, certain systems and their implementation can make using one over the other more appropriate.


Sample Rate Conversion between 32K and 44.1K I would like to know what are possible/typical efficient implementations (preferably in fixed-point DSP) of synchronous sample rate conversion between 32KHz and 44.1KHz (audio applications). Also, I'd appreciate if you can possibly share DSP utilization (MIPS, code and data memory).

Thanks!

Notes: Just wanted to make my question a bit more clear:

By "synchronous" I meant 32K and 44.1K sampling frequencies are on the same clock domain.
I don't think straight forward poly-phase implementation is a good idea here. It is very efficient for conversion between e.g. 32K and 48K, but upsampling and downsampling factors are somewhat large for 44.1K and 32K
I had an idea to do the sample rate conversion at two steps: first, convert 32K samples to 44K (upsample by 11 followed by downsample by 8, also suggested by Hilmar in the answers section here). Then do 44K to 44.1K conversion using some kind of fractional sample rate conversion (Farrow, all-pass filters, fractional phase interpolator, etc.). I used 7th order Farrow structure but I'm not quite satisfied with distortion and harmonic levels.
My application is real-time

The standard way of doing this is a poly phase filter. Unfortunately it's not straight forward and the devil is in the details.

You need to determine cutoff frequency, number of phases, length of filter, and (to a lesser extent) shape of filter
These are typically derived from your specific requirements: what's your allowable bandpass ripple, stop band attenuation, residual aliasing, sine wave modulation, up to what frequency, etc.
How exact does your final sample rate need to be? If you can tolerate a quarter percent of frequency error you can approximate the ratio as 8/11. If it needs to be better than that, you may just as well go with an asynchronous converter
MIPS and memory footprint depend very much different factors. Number of channels to convert,  $C$, filter length , $N$, and number of phases $P$. Storage requirements are roughly $(P+C) \cdot N$. MIPS for synchronous are $C \cdot N$ multiply/adds plus some overhead for phase indexing.
Asynchronous will be more, maybe $(C+2) \cdot N$, depending on how exactly phase interpolation is done.
It works well with fixed point


Will lossy compression prevent me from performing audio analysis? I'm starting a project where we will have to analyze bird recordings. We will try to measure the birds activity using the audio.

The recordings are really long (around 24 hours each) and recorded in IMA ADPCM WAV.
All files are about 1.5GB each.

When trying to import them into MATLAB, it simply locks up because of the size of the files (I think). So we were thinking about heavily compressing the files to 8KBPS mp3's, which would result in files of around 100MB.

Now I know mp3 compression is based on the psychoacoustics of the human hearing, so I was wondering whether using mp3 compression would throw away data dat would be needed for the analyzation.

If it does throw away data that we would need, what would be the best route to take? Are there other compression techniques that might be better suited to our needs?

Thanks in advance.

I would make a comment (not really sure about it) but I don't have enough rep, sorry (tell me if I need to delete this).

I just wanted to suggest trying to load the heavily compressed mp3s first, because I'm guessing matlab would decompress them to the same amount of samples and may lock up the same.

Also, if possible, I would process the files in chunks and not load them all at a time (I'm guessing matlab has a function to seek in the audio file without loading it all in the memory). You say 1.5GB each, so if there are more than a couple, you may actually be maxing your RAM (if matlab converts from 16bit to float or double, one file could probably max the ram).

Last thing, I'm not sure if you mean 8 kylobytes or bits per second. I'm pretty sure 8 kilobits per second is not gonna produce anything remotely usable.

Again, I would have made this a comment because I'm not really sure about it, but I can't :(


Cannot figure out the convolution operation in this problem A linear system S has the relationship $y[n]=\sum_k {x[k]*g[n-2k]}$, k ranging from $-\infty$ to $\infty$, between its input $x[n]$ and output $y[n]$, $g[n]=u[n]-u[n-4]$.

Determine $y[n]$ when $x[n]=\delta[n-1]$.

I know that the convolution operation can only be applied when both the operands have the same value of $n$ and that :

$$
y[n]=x[n]*h[n]=\sum_k {x[k]h[n-k]}
$$

I don't understand how $x[k]$ and $g[n-2k]$ can be convolved with each other, when they do not have the same inputs. How do I perform discrete convolution on $x[k]$ and $g[n-2k]$?

Hint:

For which values of $k$ is $x[k]=\delta[k-1]$ not equal to zero? How does that simplify the given sum?


Interpolation vs Interpolation Filter? Hi guys, I've been reading some papers on - how to remove ghosting artefacts from the Fourier Slice theorem applied to a 3D discrete image volume. The papers mention that in order to remove ghosting artefacts "employ a better(and typically wider) interpolation filter".

Currently, in Matlab I fft a 3D image space volume to Fourier space. I then extract a 2D plane at an arbitrary angle from this 3D frequency domain volume, making sure that my 2D plane passes through the centre of the volume. Since, this 2D plane is at an angle (not all points on the plane correspond to the uniform discrete values in the volume), I the use interpolation on the 3D volume in Fourier space, which I believe uses weighted linear interpolation (8 nearest point average).

What do these papers mean when they say - apply an interpolation filter, is this any different from what I am doing? If it is different, should I be multiplying this interpolation filter with the frequency spectrum and if so what does this filter look like?

Interpolation is the process of adding zeros in between your samples and then low-pass filtering to get rid of the aliases that result.  Linear interpolation is just one out of many different possible filters.  Linear interpolation filters can always be represented by FIR filters.  In one-dimension the interpolation by two linear interpolation filter is $[\frac{1}{2}, 1, \frac{1}{2}]$.  There are better low-pass filters.

It is not clear from your quote what they meant by a "wider" filter.  I'm guessing they mean "one with more taps", but it could also mean "one with a wider passband".


Convolution Offset and Factor in Frequency Domain See this link. The following Convolution in Spatial domain code is written in C#.

private void SafeImageConvolution(Bitmap image, ConvMatrix fmat)
{
    //Avoid division by 0
    if (fmat.Factor == 0)
        return;

    Bitmap srcImage = (Bitmap)image.Clone();

    int x, y, filterx, filtery;
    int s = fmat.Size / 2;
    int r, g, b;
    Color tempPix;

    for (y = s; y < srcImage.Height - s; y++)
    {
        for (x = s; x < srcImage.Width - s; x++)
        {
            r = g = b = 0;

            // Convolution
            for (filtery = 0; filtery < fmat.Size; filtery++)
            {
                for (filterx = 0; filterx < fmat.Size; filterx++)
                {
                    tempPix = srcImage.GetPixel(x + filterx - s, y + filtery - s);

                    r += fmat.Matrix[filtery, filterx] * tempPix.R;
                    g += fmat.Matrix[filtery, filterx] * tempPix.G;
                    b += fmat.Matrix[filtery, filterx] * tempPix.B;
                }
            }

            r = Math.Min(Math.Max((r / fmat.Factor) + fmat.Offset, 0), 255);
            g = Math.Min(Math.Max((g / fmat.Factor) + fmat.Offset, 0), 255);
            b = Math.Min(Math.Max((b / fmat.Factor) + fmat.Offset, 0), 255);

            image.SetPixel(x, y, Color.FromArgb(r, g, b));
        }
    }
}

Note the following lines:

    r = Math.Min(Math.Max((r / fmat.Factor) + fmat.Offset, 0), 255);
    g = Math.Min(Math.Max((g / fmat.Factor) + fmat.Offset, 0), 255);
    b = Math.Min(Math.Max((b / fmat.Factor) + fmat.Offset, 0), 255);

these lines are keeping the values of Red, Green, and Blue between 0 and 255. We are able change the values of Factor and Offset to control the contrast/brightness of the output image.

How can we achieve the same in frequency domain? I.e. how do we change the values of factor and offset in Frequency domain?

We know that a convolution operation in the Spatial domain is equal to the multiplication operation in Frequency domain. In frequency domain we are multiplying two Complex numbers. Which part of those complex numbers represent factor and offset?

Those lines that keep the range of the R,G,B channels within the valid [0,255] range work on individual pixels. There is a comparison between the pixel value and 0 / 255 levels. Unfortunately, this operation requires the exact pixel value to be known. However the Fourier transform mixes all the individual pixel values into all frequency bins and unless explicitly (and inefficiently) inverted back for a given pixel, this value is not accessible.

So it's not possible to limit the individual pixel values in the frequency domain.


derivative filter transfer function In many of the papers it is said that the derivative filter transfer function is given by: $$H(z) = \dfrac{1}{8T}\left(-z^{-2} - 2z^{-1} + 2z + z^{2}\right)$$ But no one gave the detailed information about it. Anyone has an idea about this?

This is not a standard 5-point derivative formula, the corresponding transfer function of which is

$$H(z)=\frac{1}{12}\left(-z^2+8z-8z^{-1}+z^{-2}\right)\tag{1}$$

The figure below shows the magnitude responses of an ideal differentiator (red), of the standard 5 point approximation given by (1) (green), and of the approximation in your question (blue):

As you can see, the non-standard approximation approximates the ideal differentiator only in a frequency region up to about $1/10$ of the Nyquist frequency (half the sampling frequency), whereas the standard approximation approximates the ideal differentiator quite well up to about a normalized frequency of $0.3$.

The differentiator in your question is a maximally flat lowpass differentiator. It avoids amplification of high-frequency noise. It achieves a maximally flat approximation of the ideal differentiator at $\omega=0$ and it has a given number of zero derivatives at $\omega=\pi$ (Nyquist), i.e. it is has a flat response at Nyquist, which causes high frequencies to be attenuated.

The design formula is derived in this paper. There is also a Matlab program available to design such maximally flat differentiators: link. Your specific filter is obtained by the command h=lowdiff(3,0);


Optical Flow: Method to estimate a very small displacement As I was searching for available resources( here, here, here and many more), I realize that many developed optical algorithm is to estimates LARGE displacement. However, there is very limited information or journals that put effort on estimate very small displacement. The small displacement I'm taking here is at sub-pixel level. The smaller the resolution the better.

So my question, you guys have any idea in mind that might inspire a new method, or suggest an existed method (that I might not be aware of/that can be improved)? Any method/suggestion/comment is very welcomed!

I was confronted to the same problem and I can recommend to use particle filtering. The idea is that instead of using a discretization of the space (which would favor only a limited range of motion detection), the optic flow is computed as a probability density function which is represented by a limited number of samples. These samples (or particles) are precise up to the precision of your computer and can thus detect very small displacments.

Some resources:

the original paper that led to this idea of applying the particle filter method to computer vision: http://www.robots.ox.ac.uk/~misard/condensation.html
a first implementation with details on the method: http://invibe.net/LaurentPerrinet/Publications/Perrinet12pred
a recent implementation : https://github.com/laurentperrinet/Khoei_2017_PLoSCB