Is there any fast algorithm to compute the most similar rectangular regions in two images? I am working on a pattern recognition problem to find two most similar rectangular regions in two given images. Specifically, I have two 2D (gray) images of the same sizes $I_A$ and $I_B$. Denote an arbitrary rectangular region in an image $I$ as $R^I(x_0,y_0,w,h)$, where $(x_0,y_0)$ denotes the upper left pixel of the rectangular of width $w$ and height $h$ and
$\forall i\in[0,w-1] \,\textrm{and} j\in [0,h-1],$,we have $ R^I(x_0,y_0,w,h)[i,j] = I[x_0+i,y_0+j]$.

My question is to find two rectangular regions of the same size, one in each image, (say they are $R^{I_A}(x^A_0,y^A_0,w,h)$ and $R^{I_B}(x^B_0,y^B_0,w,h)$) such that they are most similar in terms of the mean square error of theses two rectangular regions. Note it is possible to have $(x^A_0,y^A_0)\neq(x^B_0,y^B_0)$.

Ideally, I want to answer this question for all possible $w,h$ combinations. It is very costly to compute even for a given pair of $w$ and $h$. So far, I adopt the integral image technique. However, it still requires shifting image pixels. I wonder whether there is some better technique. Can anyone help?

As stated already, you can not do much to speed up a comparison between any two chosen patches. What you need to focus on in order to speed up the process is how to reduce the number of patch pairs you need to compare.

If the computational process is as expensive as I presume it is, in addition to already mentioned image pyramid, I might have another suggestion.

Extract local features and perform feature matching. This is the way to usually recognize similar images, but it works for finding objects in whole images, which means it should also match features in smaller, locally similar patches.

These guys: Sivic, Zisserman: Video Google: A text retrieval approach to object matching in videos had a good idea about spatial consistency. In short, they make sure that the groups of matches have a fairly similar geometrical layout. They use a fairly loose criterion but suggest a way to both loosen and strengthen the criterion. You do have to do a straightforward brute force matching between the feature descriptors, but that still seems far less computationally expensive than directly calculating similarities between patches.

This way, you will hopefully get clusters of matches. Then you could limit the search only to the patches that contain a certain minimal number of matches. In my opinion, these patches would have a substantially higher probability of being the most similar. This is a bit complex add on to make, but I think it could be worth a shot.


Quantization SNR of sine wave doesn't match 1.761 + 6.02 * Q I am trying to show with numpy that the quantization noise of a sine wave matches the SNR formula of SNR = 1.761 + 6.02 * Q.
The numpy code is simple:

import numpy as np
import matplotlib
from matplotlib import pylab, mlab, pyplot
plt = pyplot

from pylab import *
from numpy import *
from scipy import signal

def quantization_noise(quant):
    N=8192
    freq = 128
    x = np.linspace(0., 1., N)

    y1 = 0.5 * np.sin(2 * np.pi * freq * x)

    y2 = (np.floor(quant * (y1)) / quant)
    diff = y2 - y1

    freqs = fftfreq(N)
    x_mask = freqs >= 0

    Y1 = np.fft.fft(y1)
    Y2 = np.fft.fft(y2)

    Y1db = 20 * np.log10(np.abs(Y1) / N * 4)[x_mask]
    Y2db = 20 * np.log10(np.abs(Y2) / N * 4)[x_mask]

    plt.plot(freqs[x_mask], Y1db, 'bx', label = "input")
    plt.plot(freqs[x_mask], Y2db, 'r-', label = "output")
    #plt.plot(freqs[x_mask], Y1db, 'bx')
    plt.ylim([-140, 5])
    plt.xlim([0, 0.5])

    snr = np.amax(Y2db[3*int(freq):])
    print(snr)

    plt.plot([0.0, 0.5], [snr, snr], 'm-.', linewidth=1.0)
    plt.text(0.3, snr+4, "SNR=%4.1fdB" % snr)

    plt.grid(True)
    plt.legend(loc=1)

if True:
    plt.figure(figsize=(10,6))
    quantization_noise(8)

    tight_layout()
    plt.savefig("quantization_noise_8.png")

    plt.figure(figsize=(10,6))
    quantization_noise(16)

    tight_layout()
    plt.savefig("quantization_noise_16.png")

When I look at the results, I get an SNR of 27.4dB for 3 bit of quantization. The theory predicts 19.8db.

Similarly, for 4 bits of quantization, I get an SNR of 36.1dB: ~9dB more than for 3 bits, where you'd a delta of 6dB.

Ultimately, I want to show how with 16 bits A/D conversion, you'd end up with 98dB, but as the quantization level increases, the output spectrum gets closer and closer to the input spectrum, which is a continuous downward slope, which raises the question at which point something is considered noise instead of part of the signal.
I used applied a hanning window to better isolate the sidelobes of the main signal, which, for 3 bits of quantization makes the SNR go up from the earlier 27.4dB to 33.3dB:

I'm trying to figure out where my understanding is lacking.
How can I numerically show demonstrate the validity of the 1.761 + 6.02Q theory?
Tom

Some issues here:

Your SNR formula only applies to full scale sine waves, your sine wave has -6dB amplitude so your SNR will be 6 dB lower
The formula also implies rounding, not truncation, that's another 6 dB
You use a frequency that's a small integer divider of the sample rate, that means you are just repeating the same samples over and over again and don't get enough sample coverage to get a statistically meaningful result.
Your SNR analysis in the frequency domain is needlessly complicated and prone do imprecision and masking errors. Just do it directly in the time domain.

Here is how this would look in Matlab
%% quantization noise of a 16-bit sine wave
fr = 975.3; % something odd
n = 8192;
quant = 2^15;
% make the sine wave
y0 = sin(2*pi*(0:n-1)'/n*fr);
% quantize
yq = round(quant*y0)/quant;
% noise
yNoise = yq-y0;
% SNR
fprintf('SNR = %6.2fdB\n', 10*log10(mean(y0.^2)/mean(yNoise.^2)));

Technically, you would also have to work around the fact that the positive max amplitude of the sine might clip, but for large quantization this makes no meaningful difference.
EDIT
It may be good to review where the formula comes from and what it actually means. It all starts with the quantization noise. If we quantize and round than the quantitation noise is uniformly distributed between $[-0.5 \delta ,-0.5 \delta ]$, where $\delta$ is the quantization step. For truncation it would be uniformly distributed $[0,\delta ]$. For rounding, the resulting noise powers are
$$ P_{round} = \frac{\delta ^2}{12}, P_{trunc} = \frac{\delta ^2}{3}$$
If $B$ is the number of bits, then for a signed signal, we simply have $\delta = 2^{B-1}$ and hence we get the noise level for 16 bit rounding as
$$P_{16} = \frac{2^{-30}}{12} \Rightarrow -101.1 dB $$
That's true for ANY signal that has a reasonable wide distribution of samples. Since a full scale sine wave has a power of $-3dB$ the resulting SNR would be $-98.1 dB$


Introductory resources to mathematical properties of "dynamically recalculable" filters in audio/musical equalizers? Introductory resources to mathematical properties of "dynamically recalculable" filters in audio/musical equalizers?

Particularly, I want to understand

What mathematical features (e.g. monotonicity of coefficients) of filters are suited for real-time modulated filters?
How to analyze filters for the required features?

But I want to also understand things such as "filter instability" and the ways they are implemented in real-time (does one apply them in buffers or what?)?

Dynamic recalculation: That the filter's cutoff and resonance can be recalculated on the fly. Or that the parameters can be adjusted while the filter is running.

If you want to modulate in time your filter parameters (like cut frequency), you shouldn't use the usual DF1/TDF2 implementations derived from a polynomial fraction in z.
The reason is that they are originally LTI systems (linear time invariant), which is no longer the case now. The state that they have internally, whether an actual state or a function of the past inputs/outputs, is not reliable when the polynomial coefficients change.
Instead, you should use things like TPT schemas (see the Art of VA filter design at https://www.kvraudio.com/forum/viewtopic.php?f=33&t=350246) that will keep a state that is compatible with parameter changes. Basically, if you map an existing analog system, you model the state of the capacitors and coils and keep them continuous. This doesn't happen when using the z transform.


Recommendations for DSP Hardware I'm trying to figure out what a good DSP unit would be to try and learn how to do some embedded programming on? I would like something relatively inexpensive but that still has a decent amount of "kick" to it performance wise. I was looking at the TI ez-DSP but I don't know what to think of it.

I would like it to be programmable with C or C/C++

I would Prefer 32 bit or higher.(if possible)

Fixed point math.

Something in the Ghz clock speed.(once more if possible)

Support for hardware modulo would be a plus.

I know there are many other considerations to take into account and some of the things I have specified may not be possible to get all in one unit especially at a low cost.

I will be using it for Audio Processing and Audio Signal Generation So if a unit exists that is either compatible with ADC and DAC units or has them built in would be good.

I am open for any recommendations or any advice at all.

Thanks In advance, Alex

Since what interests you is the "embedded system" part, and since you have a low budget (this excludes anything that requires proprietary compilers), I'd recommend building yourself a board with an ARM MCU and a codec, like this one. There's less than $50 of parts - the processor, the codec and the bare minimum to get them to work.

I'm recommending this because you can easily grab a FOSS toolchain for arm-none-eabi, because these chips have a serial bootloader which allows them to be programmed with a $2 FTDI dongle, and because they have enough flash and RAM onboard for many applications - so you don't have to deal with the complications of an external RAM or flash chip. For a beginner, it's good to have a small self-contained, system you can entirely master!

You'll learn how to set up a toolchain for an embedded processor, how to use a JTAG or serial interface to flash the chip, how to configure the peripherals and hardware for your application, and you'll write code in a relatively basic environment. It'll be a fun experience! The STM32F4 is good enough to implement audio effects and 6-10 voices of audio synthesis (Used in: the Owl stompbox, Audio Damage and Make Noise Eurorack modules, DIY synth projects like PreenFM2, Audiothingies P6, or Sonic Potions LXR). UAD needs a lot of computing power because they emulate analog hardware at the component level. You are not going to do that. To give you an idea of what's possible, the STM32F4 can run about 250 biquad filters or 100 wavetable oscillators or a phase vocoder with 75% overlap. This is a general purpose architecture, but there are DSP-ish things like saturation or MAC instructions that make the generated code relatively compact on audio DSP applications.

What this won't teach you are the specifics of VLIW architectures, and the art of writing assembly for such architectures. For that, you'd need a proper DSP board from AD or TI, but this costs more, and once you build something with their devboards it's more expensive and complicated to replicate it on your own board. You can partly capture this experience with a dsPIC.

Another cheap toy you could play with: ADAU1701 board. It has a built-in AD/DA. Downside: it can only be programmed through AD's drag and drop tool, with building blocks like filters, dynamic processors... Elektor also has a spin FV-1 board - it's a very peculiar minimalistic chip that can be mastered in a few hours. Might be worth exploring if you like vintage reverbs.


Measuring sine amplitude with ADC My problem is similar to this one: Measure Sine Wave Amplitude from ADC Signal
I have a microcontroller that produces a sine wave by means of its DAC. I can choose the frequency anywhere between about 1 kHz to 10 kHz. The sine is 3 Vpp and has 1.5 V offset, so that it is an unipolar signal.
I use this sine wave to excite a passive network, and the network responds with a signal having a phase shift and a loss. What I want to do now is to estimate the loss in my network, so I need to measure its output voltage (which is also a sine wave with the same frequency).
Does it make sense to use something like a sliding FFT or a Goertzel to precisely determine the amplitude of the sine at the network's output, or is there maybe an even better approach?
My first thought after reading the linked question was that I could also use an FFT and determine the magnitude of the corresponding bin. I assume some CPU time could be saved if I implement it as a sliding FFT, or even use a Goertzel, since I am only interested in one frequency.

If you can be sure that your passive network is perfectly linear (does not create harmonics) you can just use the raw min/max values that you read from the ADC and compare them (correctly scaled) to the min/max values that you put into the DAC.
Some inaccuries will arise for sure, ie. because the DAC and ADC are not perfectly linear. If your passive network has low impedance, you would also have to take into account the internal resistance of the DAC.
If you need high precision, doing a FFT or Goertzel is the way to go.


How do I normalize accelerometer signals values to a relative scale? Different smartphone devices have accelerometers manufactured by different manufacturers.

Say I want to categorize vibrations recieved from various devices into one of the following scales :

LOW_MOVEMENT
MEDIUM_MOVEMENT
HIGH_MOVEMENT

How can I devise a scale such that -

LOW_MOVEMENT of one device is LOW_MOVEMENT of every other device when subjected to the same vibrations. And so on for MEDIUM_MOVEMENT and HIGH_MOVEMENT? (with a certain degree of possibility ofcourse)

The only data I know about the make of the accelerometer is accelerometer resolution and maximum range of accelerometer.

P.S. I plan on normalizing x^2 + y^2 + z^2 under root.

Assumptions

There is a no function in (all) target operating systems that gives you the maximum/minimum values of the accelerometer
You do not want to spend a lot of money and time to build a database of these values on your own

Suggested solutions

I assume there is no way of getting around a kind of calibration procedure. Here are two thoughts:

#1

The simplest one that I can think of is that you let people drop their phone from a certain height, let's say 1 meter. Have them to do this:

Hit a calibrate switch. This starts a count down that starts some kind of user feedback (e.g. a certain color on the screen, a tone being played) for 5 seconds.
Within these 5 seconds the phone is to be held still, approximately 1 meter above a soft ground that is safe to let a phone fall onto (e.g. blanket, bed, ...)
When the 5 seconds are over and the feedback stops (or changes, e.g. different color on the screen, or the tone stops playing), the phone is to be dropped on the soft ground.

Since you know approximately the height of the phone, you can select a part of the accelerometer data where the phone certainly is in free fall, and is therefor free of vibrations from the users hand or from the phone touching the soft surface. All you probably need is to get 10 centimeters of plain, free fall. Rotations are not a problem if you intend to use $\bar{a} = \sqrt{a_x^2 + a_y^2 + a_z^2}$. If needed, you could also let them repeat the process to gain a higher SNR from averaging. For your application, it is certainly sufficient to assume that the measured signal is due to the earth's surface acceleration of $g \approx 9.81\,\frac{m}{s^2}$. If you need slightly higher precision, you could use the geolocation of the device and an theoretical gravity approximation model.

#2

If the app is used in certain contexts and only those values are important, you could simply introduce a learning phase. Let the user wear the device/use it in the intended way, but only for learning. You could:

let the user rate the acceleration during the training phase and use this as a ground truth for your model
just collect the data and let the app figure out the levels by itself, assuming that after let's say... 13 training sessions all intensities were measured

Further ideas

In both cases, you could ask the user to share this data with the other users of your app: After a successful training session, ask them if the would share the phone model/accelerometer model (that you probably can retrieve from the OS) and the calibration values. This way, you could have the app look for calibration values online. If they are available, they could be downloaded and used without the inconvenience for the user of doing the calibration. Additionally, with such a database, you could process and thereby refine the calibration data.


Cooley-Tukey Implementation of FFT in Matlab For my course I need to implement a 30 point Cooley-Tukey DFT by transforming it into a 5x6 matrix.

I have tried to implement using the following Matlab code:

clc; clf; close all; clear all;
N = 30;
rows = 5;
columns=6;
data = linspace(0, 2*pi, N);
data = sin(data)+0.1*sin(pi*data);
plot(data);
count=1;
matrix=(reshape(data, 6, 5))'; % Read in row wise
test = matrix;
%Perform column-wise dft
for x=1:columns
    matrix(:,x)=dft(matrix(:,x));
end
%Perform Twiddles
for x=1:columns
    for y=1:rows
        matrix(y,x)=matrix(y,x)*exp((-2*pi*y*x/N)*1i);
    end
end
%Perform FFT on rows
for y=1:rows
    matrix(y,:)=dft(matrix(y,:));
end
out=20*log10(abs(reshape(matrix,1,N)));
figure
subplot(211);
stem(out); title('Using cooley-tukey');
subplot(212);
stem(20*log10(abs(fft(data)))); title('Using fft');

Where dft is defined as:

function output = dft(input)
 n = length(input);
 output = zeros(size(input));
 for k = 0 : n - 1  % For each output element
   s = 0;
    for t = 0 : n - 1  % For each input element
     s = s + input(t + 1) * exp(-2i * pi * t * k / n);
   end
   output(k + 1) = s;
 end
end

The instructions that were in our lecture were:

Load 2D array by rows
Transform each column
Twiddle each element
Transform each row
Unload 2D array by columns

Now it outputs some very different outputs:

Anyone know what I am doing wrong?

The problem is in the definition of the twiddle factors. Your indices go from 1 to columns and from 1 to rows, whereas the indices in the power of the exponential function should go from 0 to columns-1 and from 0 to rows-1. So with this line of code in the twiddle factor loop everything should be fine:

matrix(y,x)=matrix(y,x)*exp((-2*pi*(y-1)*(x-1)/N)*1i);


System Identification: smoothing in E-step of EM algorithm My question is, do the smoothed values affect the estimates obtained from the Maximization Step? If not then we could eliminate the smoothing. Once the estimates are obtained, they do not change apart from substituting the smoothed values $x_n^s, P_n^s$. So, why do smoothing?

The result of the $M$-step is expressed in terms of the matrices $\mathbf{P}_n^2$ and the vectors $\mathbf{x}_n^s$. These quantities are defined as expectations. Smoothing is used to compute these expectations, because in practice expectations are usually approximated by time averages (hence 'smoothing'). If you have another method for computing these expectations, the results of the $M$-step will of course be different, but the expressions for the results of the $M$-step remain valid, you just plug in different values for the respective matrices and vectors.


DTFT of any finite sequence in matlab I think freqz in a MATLAB toolbox, is the way to obtain DTFT of sequence.
freqz can calculate frequency response of:

H(z)=(Num)/(Den)

We can easily compute the z-transform of any finite sequence x(n) like this:

H(z)=x(0)z^0 + x(1)z^1 + ...

We know in above expression that the Den, is 1.

Recalling that: freeqz(num,den,n) gives the step response in n point.
By x be vector of x(n),

[x1freqz, x1freqzw]=freqz(x,1,3000,'whole');

must gives us  the DTFT.

1)Is it(above statement) correct? what happening if we shift our polynomial?? why?

The second way is to calculate DTFT formula completely, like this:

[X, W]=me_dtft(x1',pi,3000);
figure
title('my')
% plot(W/pi,20*log10(abs(X)));
plot(W/pi,abs(X))
ax = gca;
% ax.YLim = [-40 70];
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude (dB)')

function [X, w]=me_dtft(x,whalfrange, nsample)
    w= linspace(-whalfrange,whalfrange,nsample);
    t=0:1:size(x,2)-1;

    X=zeros(1,size(w,2));
    for i=1:1:size(w,2)
        X(i)=x*exp(-t*1i*w(i))';
    end

end

2)I confused, is the range of parameter t in above code, important?

3)Is this implementation correct? Why?

I think there must be dissonance, since the picture:

Telling us something wrong!
These transform are taken from pure sine wave(its code in in the picture), from right you can see the fft , freqz by the manner explained top and(left) the DTFT as explained earlier.

Edit after "Jason R"'s comment:
Ok, also I removed logarithmic scale, since it make me confuse.
After that, intuitively thy are alike, as you can see in next image, but why they are not exactly the same(refer to last image by logarithmic scale?)?

freqz:

[x1freqz, x1freqzw]=freqz(fliplr(XX'),1,3000,'whole');

figure
title('freqz')
% plot((x1freqzw/pi)-1,20*log10(abs(fftshift(x1freqz))))
plot((x1freqzw/pi)-1,abs(fftshift(x1freqz)))
ax = gca;
% ax.YLim = [-40 70];
ax.XTick = -1:.5:2;
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude')

Sine sample:

Fs=1000; Ts=1/Fs;
time=0:Ts:1;

Freqs=500;
Xs=zeros(length(Freqs),length(time));

for i=1:length(Freqs)
    Xs(i,:)= cos(2*pi*Freqs(i)*time);
end

XX=Xs;
XX=XX./ max(abs(XX));

figure;plot(time, XX); axis(([0 time(end) -1 1]));
xlabel('Time (sec)'); ylabel('Singal Amp.');
title('A sample audio signal');
sound(XX,Fs)

The linspace function as used is going from -fs/2 to +fs/2 for 3000 samples so the +fs/2 is getting counted twice. In contrast the FFT and freqz, which goes from 0 to N-1 do not duplicate the two end points (in this case they go from DC to 1 bin less than fs where fs is the sampling rate). Therefore the sample locations are not exactly the same leading to the difference observable in the two freqz methods.

Further you can get samples of the DTFT by zero padding as another option: fft(x, 3000).

Instead of linspace which will work just fine if you choose the proper start and stop, I like to do:

t = [0: length(x)-1]*1/fs


Frequency Axis of Discrete Fourier Transform (DFT) with Odd Number of Data Points I am trying to understand the logic behind making a frequency axis in DFT. I am using for time based light absorbance. When we have even number of data points (N= even integer), collected over a length of time A, the lowest frequency $\omega_1$ that can be resolved is 1/A, the next step is 2/A, the $k^{th}$ step is k/A, where k =1,2,3..., N/2. The highest frequency $\omega_{max}$ = N/2A based on Shannon's or Nyquist criterion. This the maximum value on the frequency axis.

How do we scale the frequency axis of FFT when we have odd number of data points (N+1) collected over a length of time A, where N is a even number? I cannot find similar reasoning for scaling when we have odd number of data points except in the book excerpt below.

The step size of frequency in case of odd number of data points:  k/A whether the data points are even or odd.

Please see the excerpt from a book called DFT: An Ownwer's Manual by Briggs. What is the author trying to say when he says that the highest frequency $$\frac{N}{(2A)}$$ does not quite coincide with the endpoint of the frequency domain which have the values $\Omega$/2?Please note that the author is using N as an even number.

Is he indicating that there is a slight error when we have a DFT of an odd number of datapoints?

This can be answered simply by considering the definition of the $N$-point DFT:

$$X_N[n] = \sum_{k=0}^{N-1} x[k]e^{-j2\pi \frac {n}{N}k }$$

where it's easy to see that the DFT just compares your $N$-point input signal $x[k]$ to a sinusoid of frequency $\frac nN$.

Thus, the lowest frequency is always $0$, and the resolution is always $\frac{f_\text{sample}}{N}$, no matter whether $N$ is even or odd. That also means that the lowest non-DC frequency is also $\frac{f_\text{sample}}{N}$.


fft function in R vs spec function from 'seewave' package? They don't give a similar frequency spectra I understand that spec will give me the frequency and the corresponding amplitude of that component, whereas fft will compute the DFT of the signal and throw the complex numbers for each component which I have to plot as per the physical frequencies.

y is the time-series signal
Fs is the sampling freq.

N = length(y);
f1=fft(y)
fa = seq(0,N-1,by=1)*(Fs/N)
plot(fa,abs(f1),type='l',xlim = c(0,3),ylim = c(0,50))
s1=spec(y,f=1000,plot = FALSE)
s1$x<-s1$x*1000   (KHz to Hz)
plot(s1$x,s1$y,type='l',xlim = c(0,3))

Red one is from spec function and blue is from fft function.

Disclaimer: I don't know even a bit of R. But I checked the documentation of the spec function: A Hanning function is applied to the analysis window.

This may very well be the source of the difference in the two plots. The fft function does not apply any window (i.e., it's only the natural rectangular window we are seeing).


Structure of function that describes sine signal I need to create a vector of sine signal.

So I'm trying to figure out what the structure of the function that describe sine signal ?

For example the function   $\sin 2x$  meets the requirements?

If the answer is no, then what is the reason?

I think you posted similar question 3 days ago, regarding your teacher claiming that $\sin 2x$ is not a sinusoidal function. Nevertheless this function is definitely sinusoidal. Otherwise how come we can say that signals can be decomposed into sums of sinusoids (and cosinusoids)? You have plenty of orthogonal waves: $\sin x, \ \sin 2x, \ \sin 3x, \ldots$, and all of them are sinusoids.

The only thing that is coming into my mind, is that:

  Discrete-time sinusoid is periodic only if its fundamental frequency
  $f_0$ is a rational number

For a sinusoid with frequency $f_0$ to be periodic, we should have:

$$ \sin[2\pi f_0 (N+n)+\theta] = \sin[2\pi f_0 n + \theta]$$

This relations is true if and only if there exists an integer $k$ such that:

$$2\pi f_0N=2\pi k $$
or, equivalently:

$$f_0 = \dfrac{k}{N} $$

To determine the fundamental frequency of a periodic sinusoid, we express its frequency as above and cancel common factors so that $k$ and $N$ have no common divisors. Then the fundamental period of the sinusoid is equal to N.

So for example:

$f_0 = \dfrac{31}{60}$ implies that $N=60$
$f_0 = \dfrac{30}{60}$ implies that $N=2$

Thus if you define your discrete sinusoid to be:

$\sin[2 f_0 t]$, or
$\sin[\sqrt{2} f_0 t]$

then these are no longer periodic in digital domain. Why? Think of that in following way: each sample at start of the new period will be shifted slightly - it's not a rational multiple of $\pi$. Below is some plot for two signals:

$\sin [2 \pi t]$
$\sin [2 \pi \sqrt{2} t]$

Sampling frequency is $10 \ \mathtt{Hz}$ and upper time limit is $8 \ \mathtt{s} $

You can see a first sinusoid is periodic - every 10 samples (1 second) you get repeating pattern. On the contrary second one, of which fundamental frequency cannot be represented as a rational number (no way to decompose $\sqrt{2}$) periods are not the same - thus your signal has no period. Check the figure below for overlay of 11 periods:


Removal of noise from vibration recordings Problem
I am analyzing accelerometer recordings that are contaminated by 'spikes' caused by issues in the digital circuitry. The spikes are problematic, because I wish to determine the peak-amplitude from these signals using an automated procedure from a large number (hundreds, perhaps thousands) of these recordings. Hence, manually determining the amplitudes is possible, but not preferred.

Background
The acceleration recordings are oscillatory in nature and reach frequencies of ~200 Hz. I am planning on using a low-pass filter in Matlab to remove the contaminating spikes in the signal. Sample rate is around 730 Hz. An FFT-based filter may be complicated by the fact that the signal gradually increases and decreases in frequency over time.

Here is a picture of an example recording including the pesky spikes (the noise) in the first 100 ms:

Question
I am using Matlab for signal processing and I am looking for advice on which digital filter for use as a high-pass filter? For example, Butterworth and Chebyshev filters are quite commonly applied if I am not mistaken, but why? Will wavelets do me any good?

If you want to detect the peak value of the time-domain signal, it is important not to change phase of each frequency signals. Otherwise, the peak value will be affected by it and will slightly change.

I think you already have recorded signal and you can process it after the recording.
If so, it is better to use zero-phased filters, like this function in matlab.

http://uk.mathworks.com/help/signal/ref/filtfilt.html

In this case, the needed function is to reduce the spike, which contains so much frequencies.
In the help page I have shown before, there is an example of filtering this kind of spikes from the relatively low-frequency signal.
I think it's what you need.

Edit

I recommend Butterworth low pass filter for its flat gain characteristics. The gain ripple of the Chebtshev would affect badly. The phase change will be corrected by using the filter backward, which the filtfit does. As you imagine, the effect of the phase change would be very little in this case, but it cannot be guaranteed.
Just for your information, Bessel would be a candidate when you do the filtering online.

If you still doubt about the effect of the phase shift (not the time-shift) by the filter, you can imagine the time-streched pulse and impulse signal in the image and code below.
These signals are almost the same but only the phase are different.
Besides the amount of the effect, the phase shift does affect to the amplitude of the signal.

I recomend you to use zero-phase filter if you want to measure the actual amplitude precisely.

clear all;
close all;

% generate signal with the same phase
N = 2^15; % number of samples for this signal
for k = 1:N/2;
    S(k+1) = exp(0);
end
S(N/2+2: N) = conj(S(N/2: -1 :2));

signal1 = real(ifft(S));

% generate signal with different phase from above
% this signal is called time-streched pulse.
N = 2^15;
for k = 1:N/2;
    S(k+1) = exp(-1j*2*pi*N/2*(k/N)^2);
end
S(N/2+2: N) = conj(S(N/2: -1 :2));

signal2 = real(ifft(S));

% and let's see them.
subplot(211)
plot(signal1)
xlim([-1e4, length(signal1)])
xlabel('time')
ylabel('all the same phase')
subplot(212)
plot(signal2)
xlim([-1e4, length(signal1)])
xlabel('time')
ylabel('only the phases are different form above')


Beam Former linear array antenna How can I design a beam former (with weigth vector $w$) that protects a signal coming from an angle of arrival $\theta_0$ and completely suppresses signals coming from $\theta_1,\theta_2,...,\theta_k$? How should I choose the number of elements $M$?

For spatial filtering there are numerous beamforming techniques and it will be difficult to compare all here, there are numerous resource available on internet.

Just to give you insight I am giving brief description of well known beamforming techniques

Phased shift beamformer also known as conventional beamformer:
This beamforming techniques tries to add constructively the signal from the provide direction. This algorithm thus shifts the signal such that they align perfectly at the provided direction and thus sums up constructively. This algorithm has the worst angular resolution amongst different beam forming techniques as all signals from other directions are also included into the result. The resolution is in the range of the beam width of linear array.
Minimum variance distortionless response (MVDR) beamformer also known as CAPON beamformer: This beamformer is based on matching the array response matrix with the inverse of the sample correlation  matrix, thus enhancing the signal from considered direction and also suppressing the effects from the interfere signals from different direction other than the direction under consideration. The resolution of MVDR is far better than phased shift beamformer.

Among these two algorithm stated above the MVDR has benefit related to separation between two signals but on the expense of more computations.
It would be of worth to also look into the Linear constraint minimum variance beamformer (LCMV).

The implementation of MVDR and LCMV weights for beamformer is present in MATLAB MVDR weights and LCMV beamformer.

  completely suppresses signals

As far i know, there is no such algorithm which can complete suppress the signal from other directions as it depends on many factors. But all beam forming algorithms can maximize the signal from provided direction and minimize the signal from the other direction. If you want the best result I would suggest also to look into eigen space based beamformer techniques.

  How should I choose the number of elements M

The number of elements 'M' provides you degree of freedom so more 'M' would be beneficial but it limits the size of antenna. Mostly you have to decide the resolution you need for the beam former and then choose 'M' accordingly it is different for each algorithm.
Higher 'M' is beneficial but you have to draw a line somewhere.


Unilateral Laplace Transform's Differentiation Property I've read in numerous places that the unilateral laplace transform is extermely useful in solving differential equations with initial conditions based on the differentiation property of the unilateral transform:
$\mathscr{L}{f′(t)}=sF(s)−f(0_−)$
What i don't understand is why this is possible to do in Laplace but not in the Fourier Transform?
Is this related to the dacaying exponential added by the Laplace transform?
I would be grateful for an in-depth answer involving intuition.

You need to look at the derivation of that property. Integration by parts gives
$$\begin{align}\mathcal{L}\{f'(t)\}&=\int_{0^-}^{\infty}f'(t)e^{-st}dt\\&=f(t)e^{-st}\Big|_{0^-}^{\infty}+s\underbrace{\int_{0^-}^{\infty}f(t)e^{-st}dt}_{F(s)}\\&=\lim_{t\to\infty}f(t)e^{-st}-f(0^-)+sF(s)\tag{1}\end{align}$$
The first term in $(1)$ is only guaranteed to vanish for $t\to\infty$ if $\text{Re}\{s\}>\alpha$ for some value $\alpha$. So even if you define a unilateral Fourier transform, that term may not vanish in general.
However, if $\lim_{t\to\infty}f(t)=0$ holds, we could use the same property with a unilateral Fourier transform. It's just much more common to use the well-established (unilateral) Laplace transform in cases where non-zero initial conditions need to be taken into account.


Wavelets and bank filters The wavelet transform has a problem as it gives poor time resolution for low frequencies and poor frequency resolution for high frequencies according to uncertainty conditions.

This appears well while using the window notation. But when using bank filters, I can't imagine this problem. So, do bank filters solve this problem?

In addition, if there is a signal with maximum frequency equal to 1000 hz, how are low and high pass filters designed to decompose the signal according to a certain mother wavelet?

Filter banks are one way implement specific discretized versions of a given continuous wavelet transform (or a given SFTF). Their performance depends a lot on how one designs wavelet voices and their subsampling.

if you to the dyadic version (or DWT), then the highest subband spans a $[1/2,1]$ (in reduced frequencies), the next one $[1/4,1/2]$ and so on. So the  frequency resolution is usually poor on the highest subband. Alternatives consist in further dividing each subband (wavelet packets), or in already splitting the data in $M$ subbands instead of $2$ ($M$-band wavelets), or a combination of the two.

Last, the filter banks are quite generic, so nothing prevents your to use an augmented one, with both wavelet-filters and STFT-filers.


Inverse Discrete-Time Fourier Transform of $X(Ω)=jΩ$ I am trying to solve it by using the properties but I can’t seem to find the same solution as on my book.

This system is a discrete-time differentiator, and in this case it's easier if you directly solve the inverse DTFT integral:

$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{jn\omega}d\omega=\frac{j}{2\pi}\int_{-\pi}^{\pi}\omega\, e^{jn\omega}d\omega\tag{1}$$

Using integration by parts it's quite straightforward to solve $(1)$. The result should be

$$x[n]=\begin{cases}0,\quad n=0\\\displaystyle\frac{(-1)^n}{n},\quad\text{otherwise}\end{cases}$$


What are the advantages and disadvantages of PN sequence over Walsh code I would like to ask what are the disadvantages and advantages of PN sequence over Walsh code?  and also Gold code over Walsh code.

I have checked online, but I didn't get clear explanation for that.

Thank you

A Pseudo-random noise (PRN) sequence ia a closer approximation to white random noise in that its energy is spread equally over the occupied frequency band (The energy is spread as a Sinc function if reconstructed with pulses just because of the pulse shape but the underlying code as a stream of impulses has a more uniform distribution), and its auto-correlation function approximates a single impulse; providing a strong correlation when the sequence is aligned with a copy of itself in time and very low correlation when there is any offset between the sequence and a copy of itself. This property makes it ideal for time alignment such as used in acquisition to find the start of a packet, or in RADAR and sounding applications to resolve a time delay with high precision. (This post details such an example with a PRN sequence: Autocorrelation to diagnose faults )  Also due to it filling an occupied frequency band evenly, PRN sequences are also ideal for training patterns to equalize a channel since the equalizer can resolve solutions only at frequencies where a signal is present.

Walsh Codes in contrast are not spread equally over frequency (as is clear if you consider the sequence of all 1's is a Walsh Code), but are completely orthogonal when time aligned. Given Walsh Codes are always an even number of digits, when you multiply one code to another in the same set and sum the digits, it will always add to 0. PRN codes generated with linear feedback shift registers (LFSR's) are always an odd number of digits, so are not able to add to complete 0 (be completely orthogonal) and further different codes can have higher cross correlations to each other. Walsh Codes are ideal for allocating users or resources in orthogonal code space when you have tight control of the time alignment of each user or resource (such as broadcasting to multiple users from a single transmitter). This orthogonality property is disrupted when the Walsh Codes have a time delay between them. This post demonstrates the channel or resource allocation and also shows the similarity to the DFT which is also simply another set of orthogonal codes: How CDMA receiver extract it's corresponding data from the receiving modulated & superposition-ed signal?

Gold Codes are generated by adding two LFSR outputs together, each generated with a separate polynomial. The advantage of this is we get many more usable codes of a given order, and the disadvantage is higher cross-correlation and sidelobes of the auto-correlation. This is what is used to generate the codes for GPS satellites and I explain Gold Codes in more detail specific to the GPS implementation at this link:

GPS Coarse Acquisition PRN Codes


audio spectrogram normalization I am performing language classification from audio signals using mel-spectrograms as my inputs for a ResNet. It works well as long as all of my audio data from different languages is from the same dataset (therefore in the same format). When I train on language A from dataset 1, and validate on language A from dataset 2, the accuracy for that language drops significantly, and instead the network guesses that a lot of the signals from language A are a different language from dataset 2.

This leads me to believe the network is first picking up on features that determine which dataset the audio came from, and then classifying the language as a secondary feature.

I tried reformatting data from each dataset so that it would contain the same amount of information. Dataset 1 (voxforge) consisted of wavfiles and dataset 2 (mozilla common voice) consisted of mp3 data, so I converted dataset 1 to mp3 data that had the same sampling rate and bits per sample as dataset 2.

I’m thinking there must be some other encoding artifact from the datasets that is throwing my network off. Has anyone else had issues with “normalizing” data across datasets or anything about converting encodings so that the information encoded is of equivalent quality?

First of all, yes, what can be achieved by any machine learning tool is only as good as the dataset it was trained on.

The ML can't guess what properties of the high-dimensional data you want it to respect and which not.

So, you'll have to make sure all your data sources are either equal in their probability distributions when it comes to the feature you're interested in (language), essentially.

A relatively easy way of doing so is random sampling, where you categorize all data sources according to languages, and then for each language, pick a subset of recordings so that all sources are equally frequent.  You then train with different randomly picked subsets.

Really, this is Bayesian statistics basics; I don't think signal processing will help you here. You can sure still try to preprocess the audio data so that it only consists of relatively recording-method-invariate features (i.e. classical audio features, e.g. Cepstrum coefficients), but then you'd basically be doing what you'd instead wanted to use a neural net for: finding appropriate properties of the signal from which to tell the languages apart, and since the problem would persist, I'm pretty optimistic that the Resnet approach would still large correlate its classification to the data source, since some information would certainly persist through that feature-finding step, and your datasets still aren't even in terms of classification target pdf.


Reversing the Order of Operators for Edge Detection? Usually, for edge detection, we perform smoothing and then pass it through difference filter. What if application of difference filter happens first and then smoothing. How the math behind the same will change?

In the classic framework both the Smoothing and the Difference Filter are applied using Convolution.
Since it is done using convolution it implies the operation is Linear Spatially Invariant (LSI).
LSI operators can be applied in any order and the result will be the same. This is also a result of the commutativity property of the convolution operator.
Let's say $ f $ is the smoothing operator and $ g $ is the difference operator then:
$$ y = f \ast g \ast x = g \ast f \ast x $$
In practice, a more computational efficient way is define $ h = f \ast g = g \ast f $.
See Difference of Gaussians.


Gesture recognition using only simple techniques I'm following a computer vision course and I have this exercise: write a program that, given an hand image, it can recognize if the hand is open, closed, in a punch, or holding an "ok" posture, using only the techniques provided till now (pixel 4/8 connected, connected region, contour finding, holes finding, blob property like centroid, area, perimeter, eccentricity, image moments, image transformation like invert/power/log/gamma correction/contrast stretching, histogram computation and equalization).

I have done it with some basic blob properties(closed hand has a low eccentricity, "ok" has a hole, open hand has a big difference between the area of the inscribed ellipse in the blob and the blob area itself with a low eccentricity).. It seems to work but the first image is a little bit problematic.

I think there could be something more to make a more robust algorithm. Maybe some kind of moment property? Could some blob axis/orientation/extreme points help?

PS test images:

You can try looking at Hu invariant moments. They can be constructed from basic moments, and are rotation, scale, reflection, and translation invariant.

Calculate them for a set of training contours first, and then apply them to the test contour.

There are implementations in Matlab and OpenCV, as far as I remember.


FM synthesis algorithm is not very efficient I have very little knowledge of DSP, I'm trying to implement an FM synthesis algorithm in code / software.

So far I've successfully implemented a FM algorithm using the following:

sine (hz * (integrate (s 0 t)))

where hz is the frequency, s is the 'signal' and t is the time. And integrate is the function I've implemented that does integration using the "Trapezoidal rule".

However if I then try to generate a bunch of samples (amplitude over time), it has to run the integrate function n times for each sample - which adds up to a lot of operations I assume and hence is rather slow.

Are there efficient algorithms I could look into? Or alternative implementations worth considering?

One approach I've considered is instead of integrating the modulating frequency on each sample, is being able to calculate these discrete and allow them to be reused between samples.

So instead of integrating:

sample 1: 0s -> 1s
sample 2: 0s -> 2s
sample 3: 0s -> 3s

Do:

sample 1: 0s -> 1s
sample 2: integrated(sample1) + (1s -> 2s)
sample 2: integrated(sample2) + (2s -> 3s)

I think you're overthinking the integrator:

Since your $s$ is a discrete signal in itself (you're doing DSP!!), the integral $S(T) = \int\limits_0^T s(t) \,\mathrm dt$ of it at any time $T$ actually collapses to the running sum!

So, as you've noticed, when you've got the value for any time $T$, the value at time $T+T_s$ (one sample later) is simply $S(T+T_s) = S(T) + s(T+T_s)$.

This should really be the least complex step in your FM.

Now, when you think about real-world FM radios, you'll find that the simple circuits (but also the more complex ones) don't actually contain an actual integrator (which you could build with an Operational Amplifier and a capacitor, for example).

What they really contain, and what makes the spectral side of things much nicer to understand, is a low-pass filter (you'll notice an integrator acts as low-pass filter). You'll notice there will be other filters in real-world broadcasting systems, so-called preemphasis filters. I'd call them a historical artifact from when receivers couldn't be trusted to equalize overly well, and that leads to a noise shaping that leads to decreased SNR (much noise) at higher audio frequencies, if you don't pre-amplify these selectively.

Now, the classical FM transmitter thus consists of an integrator-alike low-pass filter, and a preemphasis filter with a  slight high-pass characteristic. You could probably combine these two.

My conclusion to that topic always is that FM is super easy to explain superficially, but looking deep, it gets really ugly, and the way the sausageaudio is actually produced at the receiver is loosely defined by the imperfections of 70 years old radio electronics and the technological abilities of that time. On an analytical level, FM gets especially cumbersome to deal with; try yourself on the derivation of the Bessel-shaped spectra of a single audio tone in FM.  This, in fact, lead to us (as in: the communications theory basic course at university) dropping FM from the curriculum alltogether: if you'd build an audio transceiver system in this time and age, you can always get better and more resilient than FM utilizing the same power·bandwidth product by going digital. Thus, we try to avoid weighing down students with hard math and questionable applicability, even if it means that they aren't taught how the most prolific broadcasting standard of the world works.

This really shouldn't stop you, though! This is a worthwhile experiment, and building an FM transmitter is actually pretty satisfying soon as it works. Note that real-world receivers are pretty resilient to choosing wrong filters; we've actually had a bug in GNU Radio for what I guess was more than a decade in which the filter coefficients for the de-emphasis (ie. receiver side) filter were plain wrong, and no-one noticed. People simply don't expect exactly great-sounding audio over FM.


Is this statement correct from DSP aspect? So, In my report I am comparing the related works, In most of the previous work, researchers have used spectrogram as input to vanilla CNN and classify. Where I am using different handcrafted features from the signal and then feeding it to DNN.
I want to explain why the second method works better than first, for that I am comparing like this

Motivated by prior work, Several AI-based approaches have been
proposed to classify sounds using signal analysis. Most of them use
short-term magnitude spectrograms transformed from sound data as
input to Convolutional neural network (CNN). The drawback of this
method is, Since CNN is based only on a spectrogram input, some
domain-specific important characteristics might get missed in the
feature space. Such as Zero crossing rate(ZCR), Skewness, Shannon
Entropy, Fundamental frequency (F0), Formant frequencies etc

I just wanted to know, Is this statement correct from DSP aspect? because I am saying in my statement that spectrogram doesn't contain zero-crossing rate(ZCR), Skewness, Shannon Entropy, Fundamental frequency (F0), Formant frequencies etc

It's almost a matter of philosophy, i.e., difficult to argue hard facts.
On the one hand all the features you mention can be extracted from the raw signals. So in theory the network should be able to learn how to do that if they provide meaningful information for the task at hand. This is what part of the ML community is claiming: feature engineering is dead, let the networks figure out themselves what are the best features, we're only biasing them by what we think are good features, which may be wrong.
On the other hand, if you do have very good reasons to believe that some features are particularly useful, it can actually help the network. It wouldn't need to learn to discover such features first and thus, it could train much faster. This is related to the concept of known operator learning, where you embed fixed operators (that you know in advance should be beneficial to your task) into the network to make it train faster. I'm leaning towards this approach since it often makes sense.
If people used complex-valued STFTs to the input, no data were lost and all the features you mentioned (ZCR, F0, ...) could still be extracted from it. In a spectrogram, the phase is discarded, so one can argue that in fact some information is lost. If you just use your features as input though, it still does not represent the entire signal so some information is still lost and you may have to argue why you think this information is irrelevant.
Here is where it is difficult to argue hard facts and the only thing you may be able to do is show actual training results to make your point.


Detection of similar segments of two signal I have two signals, the first is a reference signal from which I selected a segment.

I'm trying to find which segment of the other signal is the most similar to this selected segment.

I'd be glad to have some ideas.

You could use correlation between the two signals or a matched filter in which you set the filter coefficients as the samples from the selected segment inverted in time. In matlab it would be something like

coef = knownSegment (end:-1:1,:);

minSize = length(other_signal) + length(knownSegment) - 1;
N = 2^nextpow2(minSize);
OUT = ifft(fft(other_signal,N).*fft(coef,N));

Then you find the maximum value of variable OUT and at that point in time the two signals are most similar


Calculating sum of a discrete signal I have to find a convolution of two signals

$h[n] = 0.5^nu[n]$

$x[n] = u[n]-u[n-3]$

the final sum, which is correct is:

$$\sum_{m=n-2}^n 0.5^mu[m] $$

note that i replaced n-k with m, that is
$ m = n-k $

So, in regards to parameter n, i have to decypher the formula.
The result is:

$$0, n<0$$
$$2(1-0.5^n*0.5), 0\le n \lt 2$$
$$1.75*0.5^n/0.5^2, n \ge 2$$

But, i don't know how to get those results.
Can you please explain.

The results are correct (maybe i just typed something wrong, if that's the case i'm sorry in advance)

One way to see it is to realize that $x[n]$ can be rewritten as the sum of 3 impulses:

$$x[n]=\delta[n]+\delta[n-1]+\delta[n-2]$$

Since $h[n]*\delta[n-k]=h[n-k]$ for arbitrary integer $k$ you get

$$h[n]*x[n]=h[n]+h[n-1]+h[n-2]=0.5^nu[n]+0.5^{n-1}u[n-1]+0.5^{n-2}u[n-2]$$

Evaluating this expression yields

$$\begin{array}{lr}0&\quad n<0\\
0.5^0=1&\quad n=0\\
0.5^1+0.5^0=1.5&\quad n=1\\
0.5^n+0.5^{n-1}+0.5^{n-2}= &\\
=0.5^{n-2}(0.5^2+0.5+1)=&\\
=1.75\cdot0.5^{n-2}&n\ge 2
\end{array}$$

which is equivalent to the result you've given.


Easiest Beat Tracking Algorithim I want a relatively easy to implement but relatively reliable algorithm for detecting beat locations in an audio file.What I am aiming to do is get the location of the closest beat to a particular time instant i.e if given current playtime, how can I get the location of next beat.Does anyone know of any algorithms that could help me solve this problem?

The matlab codes that implemented D. Ellis's algorithm are on their website:

http://labrosa.ee.columbia.edu/projects/beattrack/


How to make DFT in matlab without FFT Does anyone have a algorithm for DFT without FFT function in matlab?

and second question is how can I make zero padding in this algorithm?

As mentioned by the applesoup, you should try dftmtx(). However, if you want to write a code for generating the DFT matrix, here it is,

funtion dftmatrix = myDFTmtx(N)

    dftmatrix = [];
    for k = 0:N-1
        row = [];
        for n = 0:N-1
            row = [row exp(-j*2*pi*k*n/N)];
        end
        dftmatrix = [dftmatrix; row];
    end
end

The matrix returned by this function when multiplied with the time domain sequence or column vector, will return its DFT coefficients.


Show others how I hear myself Sooo .. I've been thinking about this stuff. We all know that we sound different from what we hear of our own voice. It is easy to find out how others hear us by recording oneself and listen to it.

But what about the other way around?

Is there a way to transform our voice in a way that others can hear us as we perceive our own voice? I find it to be a quite interesting question. Sadly I couldn't find anything on the web after a couple google searches. Has nobody thought about this or is it impossible because of some reason that I'm not seeing?

Any leads on this would be appreciated :).

The most practical attempt that I am aware of is by Won and Berger (2005). They simultaneously recorded vocalizations at the mouth with a microphone and on the skull with a homemade vibrometer. They then estimated the relevant transfer functions with linear predictive coding and cepstral smoothing.


Long burst of high amplitude noise removal I have a setup with an embedded DSPIC. This is used to link multiple audio inputs to the outputs and changed between them either when told to, or when one is active. This changing takes place on a hardware that is purchased and from what I can tell uses a relay. I have DSP software using a heterodyne filter to remove static noise as the system runs.

However, when the output switches from silent to one output on, there is a large spike in amplitude that lasts for about 3 to 4 milliseconds. This is similar to the click a walkie talkie makes.

I was wondering if there was any way to determine when this happened and then do some sort of amplitude reduction. My sample rate is only 8khz and the signal has been checked on the scope to run at 1khz during this period. I tried doing a 5 point median filter which kind of worked, but on a larger speaker system would still be loud. I also looked into doing a differentiation filter and seeing if the derivative was large. If it was I would cut the next sample point in half and continue. This had no apparent effect.

Any ideas?

If I understand you correctly, then your issue is connected with transient that is present when switching on your device, because it is a RLC circuit. Most likely these spikes are not changing over time (because values of passive components are constant), if so, then you can:

Create the template of such spike and use the normalised correlation for detection. If it's above a given threshold, then substitute samples with linear/spline or whatever interpolation.
Use the Matched Filtering approach (which is almost the same as above approach, but offers better SNR performance) and remove detected samples accordingly

If your spikes are of random nature, then there is many methods that are allowing to remove any crackles and clicks. You might refer to the literature for such solutions:

Oudre L. - Automatic detection and removal of impulsive noise in audio signals
Knapp M., Bashir R. - Audio Click Removal using LPC
Godsil S., Rayner P. - Digital Audio Restoration
Alvarez E., et al. - Detection of clicks using sinusoidal modelling for the confirmation of the clicks


Which unit does the power spectrum of microphone output have? If I sample sound with a microphone and find the absolute square of the (non-normalized) DFT, I get the discrete power spectrum (correct me if I'm wrong). Which unit does it have?

This depends a little on how thorough you want to be. Your microphone measures pressure and not actual power. That means your time domain signal has the units of Pascal and the spectrum has Pascal squared per Hz. In order to turn that into Intensity you would have to normalize by the free field impedance of air which is density*speed_of_sound and roughly 410 Newton seconds per cubic meter (at sea level & room temperature). Then the units would be Watt per square meter per Hz. You would only get to Watts per Hz if you integrate over a full sphere . This requires either an omnidirectional source or a rather tedious measurement involving a large number of microphone positions.

There is added complexity if you want to calibrate in dBSPL (sound pressure level). The 0 dB reference for this is 20 micro Pascal (which is close to the human threshold of hearing) and a pressure of 1 Pascal corresponds to roughly 94 dBSPL.


How do you properly organize data to compute multiple (independent) recursive filters at the same time taking advantage of SIMD instructions? I'm processing multiple (independent) Exponential Moving Average 1-Pole filters on different parameters I have within my Audio application, with the intent of smooth each param value at audio rate:

for (int i = 0; i < mParams.GetSize(); i++) {
    mParams.Get(i)->SmoothBlock(blockSize);
}

...

inline void SmoothBlock(int blockSize) {
    double inputA0 = mValue * a0;

    for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex++) {
        mSmoothedValues[sampleIndex] = z1 = inputA0 + z1 * b1;
    }
}

I'd like to take advantage of CPU SIMD instructions, processing them in parallel, but I'm not really sure how I can achieve this.

Any tips?

usually, SIMD instructions work by loading the values you want to apply the same operation to into a wide register, and then applying the vectorized instruction to it.

Now, collating values from all over memory to end up in one say 256-bit wide register, just to do 8 multiplications at once, probably doesn't pay.

If you, however, are able to decide the in-memory layout of these parallel streams $x_i$ to be interleaved, i.e.

\begin{matrix}
x_0[0] & x_1[0] & x_2[0] &\cdots&x_7[0]&x_0[1] & x_1[1]&\cdots &x_7[1] & x_8[2]&\cdots
\end{matrix}

then you could load the consecutive values into one SIMD register without need for any memory or register shuffling; that would be desirable.

Again, in your specific case I wouldn't expect too much from optimization for SIMD: If a single-tap IIR is already the bottleneck, your application is not CPU, but memory bandwidth bound.

If you just think you can speed up your application significantly by speeding up this especially cheap operation: I'm afraid that's rather unlikely; if it's trivial to parallelize on an instruction level, your compiler stands a good chance of doing exactly that, and really, one multiply-accumulate per output sample is pretty idyllic, so I doubt there's much to be done, considering the cost of rearranging numbers to make that work (your API defines samples to not be in interleaved order, and your compiler can't do anything about that).

There's other aspects to be looked at; among these are questions like whether double is really necessary for what is presumably fixed-point numbers with less than 32 bit of dynamic range.


Reflection Rotation matrix I'm doing a rigid body between point clouds using SVD.
Sometimes the routine produces an incorrect rotation matrix with a Determinant=-1, ie a reflection matrix.

Any idea why ?
Is there a valid way of converting the result into a correct matrix?

Take a look on Generalized/constrained Procrustes Problems.

It should be sufficient to update singular values in $\Sigma$.

The $\det(SVD)=-1$ case is discussed in more detail in 1987 paper of Arun, Huang and Blostein: "Least-squares fitting of two 3-D point sets".


Detection of signal having gaussian noise with non zero mean and a covariance matrix using Neyman-Pearson Lemma(LRT) Assume that we have a signal to be detected using the binary testing hypothesis. The received signal is distributed as follows.
\begin{align}
\mathcal H_0:& \ x = w\\[2ex]
\mathcal H_1:& \ x = s + w
\end{align}
where $w\sim \mathcal N(\mu_w,C_w)$ and $s\sim \mathcal N(\mu_s,C_s)$

I was trying to apply the NP lemma, but I'm getting stuck somewhere.

You're looking at a Generalized Gaussian Detector$^{\ \text [1]}$. The Neyman-Pearson detector decides $\mathcal H_1$ if:
$$
\frac{p(\mathbf x; \mathcal H_1)}{p(\mathbf x; \mathcal H_0)}> \gamma\tag{1}
$$
With $\mathbf w$ and $\mathbf x$ independent; you get here:
$$\mathbf x\sim
\begin{cases}
\mathcal N\left(\mu_w, \mathbf C_w\right) &\text{under} \mathcal \quad \mathcal H_0\\[2ex]
\mathcal N\left(\mu_s + \mu_w, \mathbf C_s + \mathbf C_w\right) &\text{under} \mathcal \quad\mathcal H_1
\end{cases}$$

The NP inequality in equation $(1)$ gives:

$$
\frac{\displaystyle\frac{1}{(2\pi)^{N/2}\mathrm{det}^{1/2}\left(\mathbf C_s + C_w\right)}\exp\left[-\frac{1}{2}\left(\mathbf x-\mu_s - \mu_w\right)^T\left(\mathbf C_s + \mathbf C_w\right)^{-1}\left(\mathbf x-\mu_s - \mu_w\right)\right]}{\displaystyle\frac{1}{(2\pi)^{N/2}\mathrm{det}^{1/2}\left(\mathbf C_w\right)}\exp\left[-\frac{1}{2}\left(\mathbf x-\mu_w\right)^T\mathbf C_w^{-1}\left(\mathbf x-\mu_w\right)\right]}>\gamma
$$
Taking the logarithm on both sides gives the needed test statistic $T(\mathbf x)$:
$$
\underbrace{\ln\left(\frac{\displaystyle\frac{1}{(2\pi)^{N/2}\mathrm{det}^{1/2}\left(\mathbf C_s + C_w\right)}\exp\left[-\frac{1}{2}\left(\mathbf x-\mu_s - \mu_w\right)^T\left(\mathbf C_s + \mathbf C_w\right)^{-1}\left(\mathbf x-\mu_s - \mu_w\right)\right]}{\displaystyle\frac{1}{(2\pi)^{N/2}\mathrm{det}^{1/2}\left(\mathbf C_w\right)}\exp\left[-\frac{1}{2}\left(\mathbf x-\mu_w\right)^T\mathbf C_w^{-1}\left(\mathbf x-\mu_w\right)\right]}\right)}_{T(\mathbf x)}>\ln(\gamma)
$$
You can continue from there. Also, have a look at this related question on NP detector.

$^{\ \text [1]}$: Kay, S. Fundamentals of Statistical Signal Processing: Detection Theory.


DCTs (they exist under different flavors) map data along some cosine-related frequency axis, from $0$ to... the maximum frequency of your data, which we don't know, since we don't have its sampling frequency. So the realtive frequency index is implicit, between $0$ and $255/256$ in relative frequency. $2.669825 $, mostly probably, is the DC component, here the average of your signal, divided by $\sqrt{N}$, where $N$ is your number of samples ($256$). If I try to reconstruct your data with a type-II DCT, here is what I get:

The so-called "Frequency index" should be taken with caution, as it really depends on the DCT flavor, as commented be RBJ.

I'm not sure the signal is exactly yours, as I did a guess on the DCT type. To filter above $10$ Hz, you can see to zero the amplitude of all frequency indices above some ratio related to your sampling frequency. DCTs are discrete transforms, acting on signals with "unit" sampling period. So without data sampling information, you cannot get Hertz from the DCT spectrum alone, unless you know spectral characteristics on your signal that can be detected in the spectrum.

Negative values simply mean that the cosine contribution at that frequency is negated.

But if you want to filter data, using a discrete Fourier transform is probably more appropriate.


Wavelet Transfrom + Power Spectral Density (using Matlab) I don't have background knowledge about signal processing before and new at Matlab too.

I have EEG data (with noise removed) 1x128; sampling rate = 128 Hz, It's means that I have 1 sec. data right?

OK, I need to extracted to 5 frequency bands that are Delta (0-4 Hz), Theta (4-8 Hz), Alpha (8-16 Hz), Beta (16- 32 Hz) and Gamma (32-64 Hz) by Wavelet Transform (using Matlab), What is/are function--commands should I use?

And then, I need to compute the power spectrum (PSD) for each band, What is/are function--commands should I use?

Example:

Input: X = [1:128] and I need
Output: Y = [psdOfDelta, psdOfTheta, psdOfAlpha, psdOfBeta, psdOfGamma]

Thanks.

You need to do 1-Dimensional wavelet package decomposition. Here is the matlab documentation for the command to do this.

Here is some example code based on what you've described above:

wpt = wpdec(X, 2, 'db1', 'shannon')
plot(wpt)

This will provide a nice clickable graph. As you'll see in the documentation, you can provide a lot of different parameters for the decomposition.


How to determine Nyquist Rate I have the Fourier Transform of a signal. How can I determine the Nyquist Rate of that signal? What is the formula for that?

For finite length or time limited signals, the Nyquist rate only makes sense given some finite noise floor requirement and finite total signal energy.  You look for the frequency at or above which nothing is above that noise floor and the total remaining energy is below that noise floor.  Twice that is the minimum sampling frequency to avoid aliasing causing any noise above that given floor.

For infinite length signals, there is the possibility of a frequency domain window of finite support width.


FIR filter design and implementation meaning What is the output of the convolution meaning which is used in the FIR filter?

how can conclude the output of the FIR filter in the output according to the input $x[n]$?

I don't fully understand how you are expressing your question in English, so please accept my attempt.

Discrete convolution is
$$ y[n]=\sum_{k=-\infty}^{\infty} x[n-k] h[k]
$$
where $y[n]$ is the output at time index $n$ and $h[k]$ is the impulse response.  The index $k$ ranges over $-\infty$ to $\infty$ if $h[n]$ is an infinite impulse response.

A causal convolution, i.e. the filter can only have output $y[n]$ for $x[n]$ ,$x[n-1]$ ,$x[n-2]$ $\dots$ $x[n-\infty]$
$$ y[n]=\sum_{k=0}^{\infty} x[n-k] h[k]
$$
if $h[0]\ne 0$, $h[1]\ne 0$, $h[2] \ne 0$, $\dots$ $h[N-1]\ne 0$ , $h[N]=0$,  $h[N+1]=0$, $\dots$,  $h[\infty]=0$, is of finite length, it is a FIR filter and

  what is the output of the convolution meaning which is used in the FIR
  filter ?

$$ y[n]=\sum_{k=0}^{N-1} x[n-k] h[k]
$$

  how can conclude the output of the FIR filter in the output according
  to the input x[n]?

$y[n]$ depends on a finite sum of $x[n]$  .  A FIR filter also need not necessarily be causal, but it doesn't change that $y[n]$ only depends on a finite sum of $x[n]$

I don't understand how design has anything to do with your question.


How do you model a rolling shutter in the frequency domain Say you have an image in time domain:
$$a(x, y, t)$$

The Fourier transform is then:
$$A(k_x, k_y, f)$$

So far so good. What if we now introduce a rolling shutter. How would one model that. I'd expect the rolling shutter to have an effect similar to a rect window, albeit a 3 dimensional one. Any references on this? It would be nice if I can avoid deriving it all.

A rolling shutter introduces a displacement of the lines of an image that is proportional to the relative velocity of the camera (or target) and the velocity of the shutter. In terms of "simulating" this through a video, it translates to "borrowing" a few lines of pixels from each frame to compose one final frame. Notice here how the shearing produced by the rolling shutter, which could be seen as a displacement of $x$ in 2 directions just got translated to a simple "vertical" motion in 3 dimensions.

So, all that we have to do is look at what exactly this operation does in the time domain and see how (if possible) we can replicate it in the frequency domain.

From your question, I suppose that you want to create one image (let's call it $I$) from the "integration" of the frames in $a$. Before we construct $I$, we have to make some assumptions about how $a$ was acquired. So, let's assume that $a$ was acquired with a "perfect" camera. One that does not suffer from rolling shutter itself and one that shoots at a very high frame rate. So the period of sampling $T_s = t_k - t_{k-1}$ of $a$ is $T_s \ll T_e$, where $T_e$ is the frame exposure time (the time it takes for the shutter to roll and create one image). Furthermore, let's assume that the "width" of the shutter, in terms of image lines, is $W_s$.

Then, just considering the shutter, this is indeed a "white" rectangle that moves parallel to axis $y$ across all time instances of $a$ and it looks like:

$$S(x,y,t) = rect\left(\frac{x-\frac{N_x}{2}}{\frac{N_x}{2}},\frac{y-t \times T_e}{\frac{W_s}{2}}\right)$$

This is a two dimensional rectangle. It is centred in the middle of the $X$ dimension of the frame and is as wide as the frame and in the $Y$ direction it is $W_s$ lines tall, moving downwards.

The time instances that make up $I$ are:

$$I(x,y,t) = a(x,y,t) \times S(x,y,t)$$

But if you "replay" this, all that you see is a band of pixels from each frame moving downwards.

To get rid of $t$ in $I$, you need to integrate this signal with respect to $t$. In other words, you have to sum all pixel time courses to generate just one image that is independent of $t$.

Now, this is a linear time invariant system of operations in three dimensions. To find its coefficients, you need to first transform $S_{FT} = \mathscr{F}(S(x,y,t))$.

To produce $I$ we multiplied $a$ with $S$. Multiplication in the time domain is convolution in the frequency domain. Therefore to express $I(x,y,t)$ (notice the $t$ still there):

$$I_{FT}(k,m,f) = \mathscr{F}(a(x,y,t)) \ast \mathscr{F}(S(x,y,t))$$

Here of course, notice the change of variables, $m,n,f$ represent their frequency domain $x,y,t$ counterparts.

Right, so now, we want to integrate the three dimensional complex $I_{FT}(k,m,f)$ with respect to time. We look up the properties of the Fourier Transform and realise that integration in time is equivalent to multiplication in frequency by a $\frac{1}{j \omega}$ factor plus a static DC component scaled by $\pi \delta(\omega)$. So, all we have to do now is multiply our $I_{FT}$ coefficients and offset them by these factors, in the direction of the time vector. Or in other words, in the dimension that corresponds to frames.

$$I_{FT}(k,m,f) = I_{FT}(k,m,f) \times \frac{1}{e^{2 \pi i f/N_f}} + \pi I_{FT}(0,0,0)$$

(And of course then, above, obtain the inverse to get back into "time" domain)

Of course, you can express the two as one operation and do the transformation in one step.

Hope this helps.


Detecting the rotation of digits in an image I am trying to find the rotation of a given set of digits from the image. For example:

There is no additional background information, only the digits are given in the input.

I have tried using methods as described: 1 , which uses the Tesseract library to solve for the rotation, but it does not work well in this case as there are no multiple lines (of text) involved in this scenario.

For this purpose, it can be assumed that the decimal point and the "mm" are always present in the inputs. Are there any other methods or any additional features that I can exploit to solve this problem?

You can use Hough Transform to find dominant lines in the image and then based on rho & theta parameters of the Hough transform, align your text.
First you need to remove unnecessary details from your image through closing as below:

binTextImage=TextImage<30;
close=imclose(binTextImage,strel('disk',20));

Then extract edges,
    edges=edge(close)>0;

Then apply hough transform (the normal form):

H=hough(edges);
imshow(H)

Then find the maximal point in the hough domain,

[X Y]=find(max(H(:))==H)

The rho and theta in the position X,Y represents parameters of the dominant line. Use the theta parameter to find out how much the image is rotated.

More details in MATLAB example in following links is given:
https://www.mathworks.com/help/vision/examples/text-rotation-correction.html


Lower digitisation rate for long distance communcation Why is it that usually lower digitisation rates are used for long distance communications? e.g. sound is highly distorted in undersea communications, is this because of the lower digitisation rate used?

why lower digitisation rate is preferable?

At long distances, your transmit signal loses a lot of power. As a result of that, the SNR at the receiver is possibly relatively low.

A low SNR means that you cannot transport many bits per second over that channel (Shannon limit).

Not having many bit/s means you can't push much data across, which means you must reduce the account of data in the digitized speech. It's that simple!


Baseband OFDM Implementation questions Signal processing is not my area of expertise so need help with implementing an OFDM system.

EDIT: The cabling and uncontrolled topology already exists and need a system that re-uses existing installation in large buildings

System is for wired communication on a shielded twisted pair. The topology is a combination of daisy-chaining, star-connections at the nodes which further go onto more daisy-chaining. Essentially as bad as power lines but largely interference free because of shielding and differential signalling.
Bandwidth is best kept as low as possible considering very long cable runs (half a kilometer) between nodes.
Need to be able to send up to 256kbps.

Have been reading up on OFDM from an implementation point of view. I now understand a a bit about the use IFFT/FFT to inject and extract data out of the orthogonal carriers. Need for guard intervals, frame synchronization with pilot symbols, channel estimation with pilot careers etc. Things are still coming together in my head but I have a very basic understanding of OFDM. I have some questions:

This being baseband is there a way of not having to use quadrature mixing to transmit the Real and Complex parts of the IFFT i.e. non-passband way of sending it?
I am trying to keep the demodulation part as simple as possible. Anyway I can get away with not using channel estimation i.e. non-coherent demodulation with DPSK or something similar?
Any pointers to fully implemented systems that I can dissect to learn more?
Considering the low data rate, I'll probably end up doing this on a low-end DSP.

EDIT: There is only one differential pair in the cable. The bus will have 12 to 15 nodes on it with any-to-any communication.

You're describing DSL.

So, go and implement DSL! The common ADSL is an OFDM system.

  This being baseband is there a way of not having to use quadrature mixing to transmit the Real and Complex parts of the IFFT i.e. non-passband way of sending it?

Yes, you can just make sure not to use one half of your carriers. That's mathematically not different from using a half as long DFT and and shifting it up by half-bandwidth.

  I am trying to keep the demodulation part as simple as possible. Anyway I can get away with not using channel estimation i.e. non-coherent demodulation with DPSK or something similar?

I don't understand that question. OFDM is used exactly because you'd need an equalizer to use a broadband channel, and the channel estimator is OFDM, in essence.

Yes, there's OFDM systems that use a differential modulation on their subcarriers (namely, DAB/DAB+). The rationale, however, is usually a tradeoff between achievable rate and easy implementation: DAB has to deal with doppler spread, anyway, so it would require a significant amount of channel sensing if done non-differentially; if you're wasting much SNR on channel estimation, you can as well spend 3dB of BER curve and do differential modulation, especially if you're in the early 1990's and everything else would be computationally hard.

If you're doing fixed line communication, where the channel rarely changes, using a differential modulation is very inefficient compared to having a channel sensing preamble. You'd need to build a system with a higher bandwidth to achieve the same rate you could do with channel sensing, and that in a channel where you really want to avoid being forced to use higher frequencies.

  Any pointers to fully implemented systems that I can dissect to learn more?

ADSL.

  Considering the low data rate, I'll probably end up doing this on a low-end DSP.

... end up using existing DSL link technology instead of building this on your own.

I honestly don't think the modulation is the hard part here – the link and network layer becomes interesting. Twisted pair for a couple dozen meters sounds like you actually just want to use 10/100 Base-T. Modern network cards are surprisingly good at dealing with low-quality channels, and shielded twisted pair is actually a pretty nice channel (as opposed to what you claim, similarity to powerline).

If you really don't have the chance to implement half-duplex 10Base-T, you could probably simply go, physically, for simple serial RS-485; again, the problem is likely less a physical, but network layer design effort.


MATLAB filter not attenuating enough Currently, I'm working on removing 50 Hz power line interference from an ECG signal. Before trying out notch filters, I decided to try out a simple lowpass filter with a cutoff less than 50 Hz. Here's the MATLAB code I used to make the filter:

Fs=500; %Sampling rate in Hz
Ast = 120; %Stopband attenuation
Ap = 1; %Passband ripple
Fp = 30; %Passband end
Fst = 45; %Stopband beginning
Hd=fdesign.lowpass('Fp,Fst,Ap,Ast',Fp,Fst,Ap,Ast,Fs);
d=design(Hd,'butter'); % Design a butterworth filter with the given characteristics
fvtool(d);

The magnitude response of the filter is:

As you can see, the filter has around a 180 dB attenuation at 50 Hz and even more at higher frequencies.

Now, I run this filter on my data. Here's the original data, in the time and frequency domain:

And here's the data after applying the filter:

As you can see, the attenuation at 50 Hz is nowhere near the 180 dB that the filter advertised. For a better look, here's the plot of the magnitude response computed by taking the ratio of the FFTs of the signal before and after filtering:

Clearly, the attenuation is nowhere near the level it is supposed to be at. Is this the way it's supposed to be or am I doing something wrong in the implementation?

Thanks in advance for all the help!

Edit : Here's an enlarged view of a single beat of the ECG after applying the filter. Note the 50Hz noise right after the tall peak.

(Disclaimer: I don't have a Matlab license for trying your example.)

First, your ECG signal is quite short, and you don't seem to be applying any windowing before the FFT (which is equivalent to applying a rectangular window). So your final spectrum will be the circular convolution of your signal with the window transform, and hence the tails above 40Hz.

Second, you are designing for hundreds of dB attenuation. This yields filters of very high order, which are very sensitive to numerical precision errors during calculations. In real life you think in tens of dB, not hundreds.

In general, the best approach for cleaning ECG signals is:

Notch filter at 50 Hz (a band-pass of order 6 to 8 should suffice).
Low-pass filter with cutoff around 100 Hz (notice where you stop having ECG harmonics and start having only acquisition noise).

Edit for the downvoter: The filtered signal is OK, but using FFT without care gives a misleading result. At least you should apply a smooth window function to the signal (in this case, the already filtered signal) in order to get a more faithful spectrum representation. Try with a Hamming window: plot(abs(fft(y .* hamming(length(y))))), where y is the filtered signal.


Expected value of a quantized noisy signal I am given a noisy signal $y(t) = x(t) + w(t)$, where $x(t)$ is my desired signal and $w(t)$ is the noise. In my scenario, the noise is very strong, much stronger than the desired signal $x(t)$. However, I know it is zero mean. Hence, if I can oberve $y(t)$ directly, I can recover $x(t)$ by averaging since ${\mathbb E}\{y(t)\} = x(t)$.
However, in practice I need to replace the expected value by averaging over realizations and this can only be done in digital domain. Hence, I only have access to a quantized version of $y(t)$, quantized to some finite number of quantization levels $L = 2^b$. Now my question is, can I hope to still recover $x(t)$ from averaging a quantized version of $y(t)$, i.e., is ${\mathbb E}\{Q\{y(t)\}\} = x(t)$ under some conditions on the quantizer (and, maybe the noise)?
Numerical evidence seems to suggest it works, even under very coarse quantization, i.e., when one level of the quantizer is much bigger than my desired signal $x(t)$. My intuition here is that it kind of works like dithering: the small addditive $x(t)$ alters the statistics of $w(t)$ slightly such that the quantization levels shift slightly, which allows to recover $x(t)$.
My other thought in this direction is that under some conditions, quantization can be treated like additive noise and the randomness of $w(t)$ may lead to noise realizations being uncorrelated over realizations in the ensemble.
Can we back this intuition by a more rigorous analysis and what would be the right tool for this? Are there limits to when and how this can work? Are there other reasons that could prevent this from working in practice?
In case the background matters: I need to detect very weak signals in measurements but I can repeat the measurements many times and each time, the signal will appear in the same spot (in time) so in theory even if in each measurement my desired signal is much below the receiver's sensitivity, I should be able to recover it by averaging enough trials. I think people do things like this in areas where we need to go to the extremes regarding sensitivity, e.g., astronomy/physics and the like.

This will should work fine.
For each realization $k$, we can write
$$y_k[n] = x[n] + w_k[n] + q_k[n] = x[n] + v_k[n]$$
So basically define the "effective" noise $v[n]$ as the sum of the analog noise and the quantization noise $q[n]$. Averaging will converge towards the desired signal $x[n]$ as long as two conditions are met:

$v[n]$ is uncorrelated with $x[n]$
$v_k[n]$ is uncorrelated with $v_m[n]$ for $k \neq m$

So as long as the quantization doesn't create any correlation with the desired signal or across realizations, your good to go.
In general this will always be the case if the analog noise is significantly larger than the smallest quantization step and also uncorrelated across realizations. This will even work with a 1 bit quantizer, it just going take a really long time to get any decent signal to noise ratio.
EDIT
The first condition will NOT be met if the analog noise is smaller than the quantization noise, i.e. $  \sum w_k^2[n] < \sum q_k^2[n] $
In this case the quantization noise dominates and is indeed correlated to the signal. You can just dither it: add more noise to decorrelate. The choice of the dither is primarily determined by the quantizer and not by the signal itself. If you know roughly the spectrum of your signal you can throw in a noise shaper as well to help things go faster.
It's a bit counterintuitive: by adding noise, the final effective quantization noise gets smaller.


What function to apply to a signal to expose its burstiness (i.e., make the burst part more evident than the non-bursty part) I am working on some signals one of which I have attached below:

What can I do (i.e. functions that I can apply to this signal) to get something like the signal below:

Note how one can draw a line at 0.25 to differentiate between rownumber <13000(or close), and rownumber>13000

-currently, I have tried checking how far each value of the signal is away from the mean, but that didn't seem to help either (i.e., the signal still did not look anywhere near my desired). By the way, I am "not" a signal processing guru, which means I probably would appreciate explanation intended for someone with little experience in signal processing.

Edit: Based on Dan's answer, I decided to clarify something on the question. I am using this signal for a machine learning related research, and the assumption is that I only have the sparse area, but this figure helps give insight into the nature of the dense part of the signal. Hence, in the above paragraph when I said I tried the deviation from the average, I mean as I received new signals (which looks like the dense part shown in the fig), I tried to find out how far it is from the mean of the sparse part of the signal. So I guess my new revised question will be- if you have the sparse part of the signal, and you want to use its nature, shape or amplitude to differentiate between the dense (which at the initial time is unknown) and the sparse signals, how would you do this?

You can apply a simple moving average filter and based on the result of the filter make a decision if you are in the sparse area or dense area, as your average will certainly be larger in your dense area. There is a resolution vs false alarm trade space that you will need to make; your result will be more accurate the longer you can make your moving average filter but your precision in determining where the transition occurred wil be reduced.

So in this sense, to compare to what you tried, instead of determining how far each signal is from the mean, take the mean over shorter intervals (as you move through the sequence), and based on the result of the mean compared to a threshold, decide if you are in the sparse or dense area, and if your decision is sparse, you can clip the result as you have done.

An alternate and simple algorithm would be a window counter, where you allow a certain number of counts (window) between samples above a certain threshold: if the sample comes in above threshold before the count is up you declare a positive, and if not you declare a negative and the count in reset. Further, a certain number of consecutive positives or negatives is required to change state between dense and sparse (this adds hysteresis and avoids a chatter). I don't know the detailed characteristic of your dense to suggest actual count intervals, but this is something you can consider playing with.


How to prove a filter has linear phase response? I have designed an FIR filter to have linear phase response using odd-symmetry design. The coefficients of this filter are {2,1,3,1,0,-1,-3,-1,-2}. I am now being asked to prove it has linear phase response.

I don't think saying odd symmetry design is deemed a sufficient answer.

Remember that there are four types of linear-phase FIR filters. The filter in your example is a type III filter: odd filter length and odd symmetry. The frequency response of such a system has the form

$$H(e^{j\omega})=A(\omega)je^{-j\omega(N-1)/2}\tag{1}$$

where $N$ is the filter length (number of taps), and $A(\omega)$ is real-valued and odd, i.e., $A(\omega)=-A(-\omega)$. Now it's up to you to show that your filter has the form given by $(1)$. Make use of

$$\sin(x)=\frac{e^{jx}-e^{-jx}}{2j}\tag{2}$$

Note that because of the factor $j$ in $(1)$ and because $A(\omega)$ is odd, there's a phase jump of $\pi$ at $\omega=0$, so the phase is not strictly linear. Such a phase response is called generalized linear phase.


$\tt fft$ vs. $\tt fftshift$ in MATLAB: Please explain I have been looking around on the internet for a really clear explanation for why fft and fftshift give different results.

I have a very long time series for which I have applied a Hanning window. The sample rate is $8\textrm{ Hz}$ and the recording is $\approx 9$ hours long ($262144$ point time series). The time series is assumed to be stationary.

When I do an fft in MATLAB, I get the following result:

When I do an fftshift in MATLAB, I get the following result:

Four questions:

How do I decide which one to use?
Which is "correct"
Why are they giving such different answers?
Why is the fftshift only real-valued while the fft is complex?

Pseudo-code:

t=load('t.mat'); %load time series
N = length(t);
dt = 0.125; Fs = 1/dt; df = Fs/N;
f = -Fs/2:df:Fs/2-df;

t_hann = hanning(t);
F_fft = fft(t);
F_fftshift = fftshift(t);
figure(1); plot(f,F_fft);
figure(2); plot(f,F_fftshift);

Any help is greatly appreciated.

Cheers

I believe you simply misunderstand what these functions are supposed to do. You get different results because they are not meant to do the same thing. I think it would be relevant that you read the documentation of these 2 functions.

Answers below:

  How do I decide which one to use?

fft : Use it when you need to compute the Discrete Fourier Transform of your time series
fftshift : Use it when you want to to swap the position of the negative and positive frequencies of a spectrum gotten with fft()

  Which is "correct"?

What do you consider correct ? fft will be correct if you need a frequency spectrum of your time domain signal. It doesn't really make sense to apply fftshift on a time domain series.

  Why are they giving such different answers?

Because they don't do the same thing.

  Why is the fftshift only real-valued while the fft is complex?

fft computes the discrete Fourier transform and by definition the output is complex.   fftshift doesn't compute anything except swaping the position of the samples, so if your input is real, you get real output.

Normally we use these two functions this way : spectrum = fftshift(fft(x))


Multiplication of fixed-point numbers I'm trying to understand the professor who wrote about Q2-fixed-point number 1.01 and 1.11 which should be multiplied. I don't understand anything.

    Multiply the Q2-fixed-point number 1.01 and 1.11 and give the result as a Q4-number. Answer:

1.01 * 1.10
-----------
 0.01    (2-komplement of 1.01)
   0
 1.1111
-----------
10.0011 (sign extension)
 0.0011

  The result of the multiplication is a positive number. Thus 1.01 is
  negated by adding the boolean complement to 1. The final result is
  0.011. This can be verified using decimal multiplication. The
  Q2-number 1.01 is -0.75 in decimal base 10 and the Q2-number 1.10 is
  -0.5 in decimal. Thus the result is 0.375 which is represented as the Q4-number 0.11.

I don't understand the above at all. I understand that 1.01 is -0.75 and that 1.10 is -0.5_10 in decimal but why use 1.10, where does it come from? Why use the 2-complement in the first place? What does the calculation mean? What does the zero mean? It is the most confusing "answer" I've ever seen.

There are errors in the presented answer. Here's the same but corrected and rewritten a bit:

 1.01 * 1.11
--------------
 carry: 1 1
        1.1111  (shift right 1.11 with sign extension)
        0.000
      + 0.01    (2's complement of 1.11)
--------------
      = 0.0011

Here, 2's complement means negation, which is done by flipping each bit and adding one least significant bit valued 1. In multiplication of signed binary numbers, if the first multiplicand (the one that is not shifted) is negative, the correct result is obtained by doing subtraction in place of the last addition by the second multiplicand, and that's the + 0.01 line. It's a binary trick that is equivalent to actually dragging along all the sign extensions, which leads to large values in the carry and more steps than necessary:

 1.01 * 1.11
--------------
 carry: 2 1
        1.1111  (shift right 1.11 with sign extension)
        0.000
        1.11
        1.1     (sign extension of the multiplicand 1.01 is ...111.01)
      + 1.
--------------
      = 0.0011

There is no need to have two binary decimals before . in the result.

The subtraction trick can be tested with different values and it passes all the tests:

              carry:             1               1
second multiplicand: 000 001 010 011 100 101 110 111
                     00  01  10  11  00  01  10  11
                   + 0   1   0   1   0   1   0   1
----------------------------------------------------
                   = 000 111 110 101 100 011 010 001
     2's complement: 000 111 110 101 100 011 010 001


Proof that module of FT of 2 independent signals is sum of modules I found on these posts (PSD subtraction and PSD of a sum of two stationary real signals)
what I expected : that, just like the variance of the sum of 2 independent signals is the sum of the variances of these signals, the psd of 2 independent signals is the sum of the psd's of these signals. I understand the explanation using the fact that the psd is the FT of the correlation function, but I am wondering if this can be prooved only using FT-space calculation.
In other words, how to show that, for $x$ and $y$ being 2 independents discrete real signals, then:
$$ |FT(x+y)|^2 = XX^* + YY^* + XY^* + X^*Y = |X|^2 + |Y|^2 + XY^* + X^*Y = |X|^2 + |Y|^2 $$
where $X = FT(x)$, $Y=FT(y)$, i.e how to show that :
$$XY^* + X^*Y = 0$$
(using DFT notation for eg, x and y being N-lenght signals).

After a few try, I end up here, for spectrum variable k, DFT on N points :
$$(XY^* + X^*Y)[k] = \sum_{i,j}x[i]y[j]2\cos(2\pi k/N(i-j))$$
which really is 2 times the real part of the CSD. But then I don't see how to that this equals 0 using the fact that $x$ and $y$ are independents (which I didn't used so far in my calculation).

Cheers

Assume $x[n]$ and $y[n]$ represent two jointly WSS (wide-sense stationary) discrete-time random processes. Further assume that they have zero means, without losing the generality.

If $x[n]$ and $y[n]$ are independent, then they are also uncorrelated, and since they have zero means, they are also orthogonal; i.e.

$$ E\{ x[n]y^*[n-k]\} = 0 \tag{1}$$

Let $$z[n] = x[n] + y[n] \tag{2} $$ be the sum process, then PSD (power spectral density) of $z[n]$ is

$$ S_z(\omega) = FT\{ r_z[m] \} \tag{3} $$
where  $r_z[m]$ is the Auto-Correlation sequence of $z[n]$:

$$ r_z[m] = E\{ z[n]z^*[n-m] \} \tag{4} $$

Based on Eq.1 and Eq.2, expand Eq.4 as:

$$
\begin{align}
r_z[m] &= E\{ z[n]z^*[n-m] \} \\
& = E\{ (x[n]+y[n])(x^*[n-m]+y^*[n-m]) \} \\
& = E\{ x[n]x^*[n-m]\} + E\{ y[n]y^*[n-m]\} + \\
   &+  E\{ x[n]y^*[n-m]\} + E\{ y^*[n]x[n-m]\}  \\
& = r_x[m] + r_y[m] +  0 + 0 \tag{5}
\end{align}
$$

The zeros on the last line follow because of Eq.1, ($x[n]$ and $y[n]$ are othogonal).

Hence what follows is that the PSD of $z[n]$ is the sum of PSDs of $x[n]$ and $y[n]$ :

$$ S_z(\omega) = S_x(\omega) + S_y(\omega) \tag{6}$$

Everything up to this point is theory based on mathematical (probabilty) model of random processes. What you are after, however, is about the practical computation of the PSD estimation and how the independence of $x[n]$ and $y[n]$ is reflected into tho Fourier calculations.

Assuming we use the periodogram estimate for PSD :
$$ \hat{S}_x(\omega) = \frac{1}{N} |X(\omega)|^2 \tag{7}$$

then what we want to show is, as you have also worked out,

$$ \hat{S}_z(\omega) = \hat{S}_x(\omega) + \hat{S}_y(\omega) \tag{8}$$

which implies that

$$ |FT(x+y)|^2 = (X + Y)(X^*+Y^*) = |X|^2 + |Y|^2 + XY^* + X^*Y = |X|^2 + |Y|^2 \tag{9} $$

and Eq.9 implies that:

$$ XY^* + X^*Y = 0 \tag{9}$$

Eventually it's convincing to show that $XY^* = 0$.

To show this is you can use the ergodicty assumption. If $x[n]$ and $y[n]$ are ergodic in their cross-correlation, then the theoretical expectation can be replaced by the practical sums as:

$$E\{x[n]y^*[n-m] \} = \lim_{N \to \infty} \frac{1}{N} \sum_{n=0}^{N-1} x[n]^*y[n-m] \tag{10}
$$

then it can be shown that the estimate is equivalent to convolution of the sequences :

$$
\frac{1}{N} \sum_{n=0}^{N-1} x[n]y^*[n-m]  = x[n] \star y^*[-n]\tag{11}
$$
But from Eq.1 we know that the $x[n]$ and $y[n]$ are orthogonal implying that

$$ x[n] \star y^*[-n] = 0 \tag{12}$$

finally, if this sum is zero then it's Fourier transform  also is zero, completing the proof.

$$ x[n] \star y^*[n] = 0 \implies X[k]Y^*[k] = 0 $$

Note that, as MArcus has already commented, the equality holds only for the ideal case of infinite length process realization and for any practical finite length sample sequence, the Eq.s (6,7,8,9,10,12) won't hold true;i.e., they will be approximate).