texts: Ultrasonic speaker calibration and emitting calibrated signals

texts:
I'm trying to calibrate an ultrasonic speaker with the aim of emitting predictable signals, but alas I keep running into trouble, probably due to my lack of DSP-fu.

A little background

I want to be able to playback sounds as close as possible to a calibrated recording I have. As far as I understand the theory, I need to find the speakers transfer function and deconvolve the signals I want to emit with it. Something like this (in the frequency domain):

X -> H -> XH

Where X is the emitted signal H is the speakers transfer function and XH is X times H. A division (./) should now give me H.

Now, in order to emit a calibrated signal, it should be divided by H:

X/H -> H -> X

What's been done

Placed speaker and a calibrated microphone 1 m apart on tripods.
Recorded 30+ linear sweeps 150KHz-20KHz, 20ms long, and recorded @ 500 KS/s.
Aligned and averaged signals with the Matlab/Octave script below, under the script the resulting signal can be seen.

files = dir('Mandag*');

rng = [1.5e6, 1.52e6];

[X, fs] = wavread(files(1).name, rng);
X = X(:,1);

for i=2:length(files)
    [Y, fs] = wavread(files(i).name, rng);
    sig = Y(:,1);
    [x, off] = max(xcorr(X', sig'));
    off = length(X) - off;
    if(off < 0)
        sig = [zeros(1, -off), sig(1:end+off)'];
    elseif (off > 0)
        sig = [sig(off:end)', zeros(1, off-1)];
    end
    X = X + sig';
end
X = X/length(files);

Fourier transformed X and XH and did the  calculations mentioned above, the result looks plausible. Below is a normalized plot of H (purple) and X/H (green).

The plot has been truncated to the relevant frequencies.

Please let me know if I'm going about this the wrong way.

My question

After calculating X/H I need to transform it back to the time domain, I assumed this would be a simple ifft(X./H) and wavwrite, but all my attempts so far have failed to get any plausible answer. A frequency vector Hf, H, and X can be found here in mat7-binary format.

Maybe I'm just tired and there is a simple solution here, but at the moment I can't see it. Any help/advice is much appreciated.


texts: Symbol synchronization with unknown symbol rate

texts:
I'm learning about the early late gate symbol synchronization and understand that it needs to be given the symbol rate ahead of time.
How could a multi-symbol rate system use an early late gate?
It is my understanding that if the rate is known, the error is computed by going through the signal every samples-per-symbol and computing the error and eventually the timing offset converges.
But if the system has multiple symbol rates, how can a receiver know what the symbol rate without a control channel?
The one idea I have is to first use a symbol rate estimator, but the problem is that there will be some error from the estimator, then the error gets used in the symbol synchronization. Ideally, there is some synchronization method that don't need to know the symbol rate ahead of time or one that is resistant to small errors in the symbol rate estimate.
Edit
This diagram is what I'm looking at:

The sampler block is the one which needs to know the symbol rate since it is taking samples every $T_{sym}$ seconds. But if the symbol rate estimator is off, then even if the timing offset $\delta$ is estimated correctly, the sampler will take samples like $t=\delta$, $\delta+\hat{T}_{sym}$, $\delta+2\hat{T}_{sym}$, ... . Is there a way to detect that the sampler is off and refine the symbol rate estimate?


texts: Modulation with Maximum Likelihood decoder and capacity

texts:
My question comes from the high spectrum efficiency claim of OFDM. This seems not correct IMHO, at least theoretically.

For me, we lose energy in using cyclic prefix, hence OFDM should be a waste of theoretically achievable rate. That said, OFDM has advantage in implementation and it is easier to get a high throughput.

Long story short, my question is: if complexity is not my concern, can (any) modulation (without channel coding) combined with a Maximum Likelihood (or Maximum a Posteriori) decoder achieve asymptotically the Shannon capacity ?


texts: What is the difference between Statistical Signal Processing and Estimation and Detection?

texts:
Looking at the syllabus of Statistical Signal Processing in different university I see a lot of correlation with that of Estimation and Detection? In some universities, these are seen the same. For example

Rice ELEC 531: Statistical Signal Processing
MIT 6.432: Stochastic Processes, Detection, and Estimation

Now, I'm confused if these two courses are separate courses or can be seen the same. It seems Estimation and Detection is a more specific and/or old naming for Statistical Signal Processing. Is this correct? Do you know of any university offering both of them? If so, what are the differences?


texts: What's the number of FFTs in Welch-based PSD?

texts: I've read from the book "Proakis, J. G., & Manolakis, D. G. (2004). Digital signal processing. PHI Publication: New Delhi, India." Chapter 12 that the number of FFTs required for estimating the PSD from Welch method with 50% overlap is $2N/M$ where $N$ is the amount of data that you have from the signal, $M$ the size of the segments you choose for the method and I guess the $2$ is for the 50% overlap but from the articles I've read they usually do one FFT less than the one described in the equation, why is that? or what am I missing from this interpretation? My theory is that maybe the FFT missing it's supposed to be done with the last 50% of the data and another 50% of data completed with zero padding. Am I correct?


texts: DAC output SFDR and Harmonics does not meet at particular decimal frequencies

texts:
We are using a current output DAC to perform waveform generation in the band of $20\textrm{ MHz}$ to $500\textrm{ MHz}$ and sampling at $1.2\textrm{ GHz}$. SFDR and harmonics suppression requirement for our application is $60\textrm{ dBc}$. During initial testing for the above mentioned band we are able to achieve the specifications of SFDR and harmonics suppression, until outputs were generated in following manner,

For $200.3456\textrm{ MHz}$ - SFDR and harmonic suppression is $>60\textrm{ dBc}$.
For  $200.001\textrm{ MHz}$ - suppression is below $50\textrm{ dB}$.

Could anyone differentiate between the above the outputs generated?

Also these spurious levels are high at only at specific points such as: $40.001, \ 80.001, 120.001,\ 160.001,\ 200.001,\ 240.001,\ 280.000,\ldots, 440.001,\text{and}\ 480.001 \textrm{ MHz}$.

Spurious at $200.001\textrm{ MHz}$

.
Spurious at $240.001\textrm{ MHz}$

So anyone could suggest a requirements in DDS for clean spectrum at DAC output?

DDS in implemented in a FPGA (Kintex 7).My sampling frequency is 1.2GHz and it is anything to do with decimal frequency outputs such as 200.001?


texts: Why does Matlab spectrogram look stripy when windowed input is zero-padded?

texts:
I am working on a little project, and stumpled upon something I cannot understand yet. I am learning cognitive science, and don't have DSP knowledge background, only learning from what I find on the internet. I have got an 46380x1 double array, which contains my signal. I would like to draw a nice spectrogram out of it like this:
%%signal = %%too long to put here, 46380x1 double array
winsize = ceil(length(signal)/10); %% =4639
win = rectwin(winsize);
shift = ceil(length(signal)/150); %% =310
nfft = winsize;
srate = 78;

spectrogram(signal, win, winsize-shift, nfft, srate, 'yaxis');
ylim([0.015 2]);

It works fine, and shows this:

But if I would like to zero pad the windowed signal (Matlab docs say, increasing nfft will do that), the output gets these stripes in it, and no matter what factor I scale "winsize" by, It will always show some kind of distortion in the spectrogram. For example I set nfft=round(1.5*winsize) before making the spectrogram, and the its output changed to this:

Why is this happening? I know I cannot increase the amount of information in the signal by zero padding, but it might be more readable to the human eye if not so pixelated, that's why I am trying to apply zero padding here. Can you please help, what I should try to achieve this without getting those lines?
Anyway, there is a huge magnitude output component of the spectrogram, throughout the whole time axis, at the very close-to-zero frequency, if I view the trajectory starting from 0hz. view(-135,65); ylim([0 2]); Why is that?

Thanks in advance!


texts: Phase difference measurement of a signal sampled with two different sampling frequencies

texts:
I am working on phase interferometry for locating a transmitter. The direction of arrival of an incident wave can be estimated from the phase difference caused by the antenna separation as shown

In order to compute the relative phase difference between an incident wave at both receivers, their sampling rate should be the same and both receivers should be phase-matched.
In my case, the sampling frequencies on both receivers are different. The phase difference computed in this scenario will be.

$\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad   \Delta\phi = \phi_o + \Delta \omega t $

where $\phi_o$ is the initial phase, $\Delta \omega$ is the angular frequency difference between the two sampled signals due to different sampling frequencies and $t$ is the time instant.

It is clear that the phase difference will vary with time and frequency. The sampling frequencies for both the receivers are $737MHz$ and $631MHz$. I am using 256pts complex FFT for phase computation.

Is there any single shot solution for comparing the relative phase difference between the signals received with different sampling frequencies as mentioned?


texts: EEG signal windowing and framing in python

texts:
I would like to do EEG signal(using 1 channel) Windowing and Framing.

Actually I am doing my thesis and following a research paper.In this paper they explained

"Windowing consists of dividing entire EEG signals into a number of windows of specific time intervals which is 1 min selected. Framing consists of the dividing selected window into a number of frames of the specified time interval. Frames can be overlapping and non-overlapping. Typically it is selected as 1sec."

By following this paper I have done my implementation below:

data,times = Be_band[:]
print('data shape ==>row,colm',data.shape)  #data shape ==>row,colm (1, 2561)
print('time shape ==>row,colm',times.shape) #time shape ==>row,colm (2561,)
print('data len ==>',len(data))              #data len ==> 1
fs=128  #sampling frequency

Then I perform the windowing and framing task as follows:

window_hop_length= 0.01 #10ms
overlap=int(fs*window_hop_length)
print(overlap,'Overlap') # got  1
window_size=0.60 #1min in window size,
framesize=int(window_size*fs)
length= len(data)
number_of_frames=int(length/overlap)
frames=np.ndarray((number_of_frames,framesize))

print(framesize,'frame size')
print("number_of_frames=",number_of_frames)
print( frames)

And I got the result of my implementation:

1 Overlap
76 frame size
number_of_frames= 1

The number of frames make me confused,because it is one.

Can any one help me how can I implement this correctly in Python?

For coding I used this as a reference.


texts: How to convert raw pulse sensor data into actual pulse?

texts:
I've built a simple device, which sends readings from a pulse sensor to a PC.

I get these data:

What technique or algorithm can I use in order to calculate the pulse rate (beats per minute) from these data?

Update 1 (03.06.2013 21:01 MSK): The time series represents my heart beat during approx. 3 minutes. According to Wikipedia, adults have a pulse rate between 60 and 100 beats per minute.

Update 2 (09.06.2013 23:53 MSK): Here you can find signal values, which were measured 1000 times per second during approx. 5 minutes.


texts:
One pertinent attribute is normalization, which measures performance in describing signal spectral amplitude and energy, like here. Others are robustness to noise, time vs frequency resolution. Anything else? And are there examples of each being computed (coded, discrete) on test signals?
The goal's to compare Continuous Wavelet Transform implementations, but I don't have much for comparing except norms.


texts: Getting introduced to signal processing for remote sensing

texts:
I work mainly in GIS (geographical information systems) but I would like to deepen my knowledge of remote sensing (analyzing imagery acquired remotely from airborne or satellite sensors). I was advised to learn as much about signal analysis as possible. I am assuming that signal processing and analysis refer to the same techniques.

Because signal processing seems to be applied to all sorts of different areas, I'm not quite sure where to start. I should specify that beyond two university calculus courses (on which did quite well), I do not have a background in physics or math.

Could anyone help me plan a sort of "road map" to learning signal processing, preferably applied to remote sensing? Thank you.


texts: Why non gaussian variables are independent

texts:
I have read a lot about ICA. I think I could not find the answer that why non-Gaussian Variables are independent.  What I understand, Central Limit Theorem states that

  distribution of sum of independent variable tends toward more Gaussian than its original random variables.

$
s_i
$ are the number of original independent sources in ICA; whereas ICA model is $x=As$. So we can define $y=w^Ts$. The main goal is to find the unmixing matrix $w$ that maximize the non-Gausaanity. So my question is what is non-Gaussanity here and why its necessary to maximize it to extract the original sources.

Any Enlightenment Please!


texts: scipy.signal.wiener for audio processing

texts:
Does some have an example of what a Wiener filter (that can be used directly with scipy.signal.wiener) can be useful for, in sound processing (it seems that such adaptive filters can be useful for noise reduction, etc.) ?

I tried with various (noisy + sinusoids) soundfiles (read into an array x) with the command

y = wiener(x)

But it doesn't give such great result... (Normal : I cannot expect something magical with just a simple function that doesn't take any parameter!)

I just wondered what this function could be useful for (sound processing).


texts: Calculate the Output of Linear Time Invariant System Given it Impulse Response

texts:
A filter is defined as $ h \left[ n \right] = \delta \left[ n \right] - \delta \left[ n - 1 \right] $.

Given a signal $ h \left[ n \right] $ defined as:

$$ x \left [ n \right ] = \begin{cases}
1 & \text{ if } x \geq 0 \\
0 & \text{ if } x < 0
\end{cases} $$

Let $ y \left[ n \right] = \left( x \ast h \right) \left[ n \right] $.
What is the value of $ y \left[ -1 \right], \, y \left[ 0 \right], \, y \left[ 1 \right], \, y \left[ 2 \right] $?


texts: Calculating phase shift between two sines of equal, unknown, amplitude and frequency

texts:
I've got a source signal $s_1(t)$ and a system affecting that, yielding a phase-shifted version $s_2(t)$;

\begin{align}
\color{blue}{s_1(t)} &= A\sin(t+\phi_1)\\
\color{red}{s_2(t)} &= A\sin(t+\phi_2)\\
\end{align}

The amplitude $A$ is unknown. How many samples (sampling happens periodically) would I need to verify that $\phi_1-\phi_2 = \frac32\pi$?

As a practical example, given these 2 points on each curve: on the purple curve, the original signal:

point A (231,3.533151)
point B: (230,3.614015)

and on the the filtered signal: the red curve:

point C (231,1.398873)
point D (230,1.174265)

How would you confirm the phase shift of $\frac32\pi$, given an equal -unknown- amplitude and equal -unknown- frequency of the 2 signals? Is a numerical solution possible for online, causal analysis? Or do you need 3 points on each curve? Thank you.


texts: FFT magnitude of ADC samples, converting back to $\textrm{dBm}$

texts:
I am generating a sinusoid from a signal generator set at a specific $\textrm{dBm}$ and inputting this into an ADC as a real input (not IQ). I take a 256 point FFT on the resulting samples (raw signed ADC decimal numbers, not voltages). I plot the FFT magnitudes and see that I have a single spike (well 2 spikes because of redundant/symmetry) that has a certain unsigned-decimal magnitude. The FFT results are basically the same as performing fft() in MATLAB.

I'd like to be able to take this one single data point (FFT magnitude at peak bin) and produce an estimate of the original $\textrm{dBm}$. I cannot just take the IFFT to reproduce a time domain sequence.

This is my attempted method, but the results are coming back too high:

ADC Termination = 50 Ohm
ADC reference voltage = +/- 1.25V
ADC number of Bits = 12

Normalize FFT Magnitude : mag_n = mag/256
Apply ADC scaling: volts = mag_n * (2.5/2047)
Power at Termination Load: power = (volts^2)/50
Convert to $\textrm{dBm}$ : dBm = 10*log10(1000 * power )

I would assume that if I got the math right, my results would come back too low since some of the spectrum energy and smeared into the side lobes and not accounted for with my single data point. Any hints?


texts: How is LMS/FXLMS noise cancelling different than simple polarity inversion?

texts:
Consider a noise-cancelling headphone:

If I have a noise signal from the outside world, mic it, flip the polarity by multiplying by -1 (with an op-amp or digitally), delay by the appropriate amount, and playback out of the speaker -- the summation of the noise and anti-noise at the ears should approach zero (with better performance at low frequencies).

How does FXLMS or LMS improve upon this? Why is polarity inversion ('phase cancellation') not as effective at cancelling noise, at least in the literature (e.g. Kuo 1994)?


texts:
I'm working on rolling noise emissions from cars. I have a model for rolling noise emissions which gives me sound pressure levels in third octave bands (29 in total, between 20Hz and 8kHz), depending on vehicle speed, and road/tyre combinaison.

I'd like to fit a real recording of tire/road noise to the model, while keeping the spectral properties of the recording. The signal of rolling noise is mainly stochastic

In blue, the spectrum of the recording using Welch's method, and in red, the levels i'd like to reach.

What kind of methods can i use in this particular case?

Thank you for your help!


texts: Decoupling two signals where one is a delayed version of the other

texts:
I have a signal $y(t)$ that consists of the sum of two signals where one is a time delayed version of the other.

$ y(t) = x(t) + x(t-t_{d}) $

$x(t)$ is known but the time delay $t_d$ is not.

I would be interested to know what (if any) techniques could be employed to recover the time delayed signal from $y(t)$.

I believe for techniques such as ICA it is required to have two receivers to decouple the signals fully however neither signal would be known in advance for this case.


texts: Number of microphone limitation for Spherical Harmonics based DOA estimation

texts:
I am following this book "Theory and Applications of Spherical Microphone Array Processing", as well as this paper on DOA estimation "Eigenbeam-ESPRIT for DOA-Vector Estimation".

My dataset is composed of a small number of channels (4) from the Eigenmike32. The book states "all spatial sampling schemes require at least $(L_f + 1)^2$ microphones to sample a sound field of order Lf without aliasing.", where $L_f$ is the order of the sound field (3.4 Spatial Sampling). Both the paper and the book use 32 microphones in their implementations.

First, I do not find any information on what is the order of the sound field. It sounds to me that this is a signal-dependent quantity, much like a sampling rate is lower bounded by the highest frequency of the signal you want to sample. Did I get this right? Can anyone elaborate please?

Second, does that mean that I can do very little with 4 microphones? Is there a lower limit on the number of microphones for a DOA estimation and for a number of active sources estimation?

What other, none neural network, methods can I use to estimate the number of active sources and DOA with (no necessarily by the same model)?


texts:
Any idea on separating multi-panel figures?

There are 2 different sub-figures connected each other.

Thanks!


texts: What Does It Mean by Gradient / Divergence / Curl of a 2D Image $ g \left( x, y \right) $?

texts:
Suppose, a 2D image is represented by g(x,y).

Then, what does it mean by "Gradient of g(x,y)", "Divergence of g(x,y)" and "Curl of g(x,y)"?

Please, avoid formal math for explanations.

Please, use layman's terms and informal discussions.


texts: Feature vector for music recommender system

texts: For building a music recommender system(content based) which signal features should i use as feature factor ?? Will it be peridogram,spectrogram or fft or anything else ? This system will recommend music corresponding how much they are similiar. So which features correspond strongly to how much signals are similiar ? Assume correlation measures signal similiarity (but i havenot seen any implementation based on correlation on signals).


texts: how to cut off frequencies above 10 Hz with Welch method

texts:
This is the first time I have worked with any code so I am rather unfamiliar with all the lingo and processes involved. I am trying to convert a signal to the frequency domain and want to find the spectral power within the 0.5-1.0 frequency band. I have the code to produce the relative power within this band, but I know that above 10 Hz it is largely noise and I would like to remove it so it does not interfere with the relative power calculation. This is the code I have been using.

import numpy as np
data = np.loadtxt('Post_EO_Foam_CTRL.txt')

import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font_scale=1.2)

# Define sampling frequency and time vector
sf = 200.
time = np.arange(data.size) / sf

from scipy import signal
# Define window length (4 seconds)
win = 4* sf
freqs, psd = signal.welch(data, sf, nperseg=win)

# Define prop lower and upper limits
low, high = 0.5, 1
# Find intersecting values in frequency vector
idx_prop = np.logical_and(freqs >= low, freqs <= high)

from scipy.integrate import simps
# Frequency resolution
freq_res = freqs[1] - freqs[0]  # = 1 / 4 = 0.25
# Compute the absolute power by approximating the area under the curve
prop_power = simps(psd[idx_prop], dx=freq_res)
print('Absolute prop power: %.7f mV^2' % prop_power)

# Relative delta power (expressed as a percentage of total power)
total_power = simps(psd, dx=freq_res)
prop_rel_power = prop_power / total_power
print('Relative prop power: %.3f' % prop_rel_power)


texts: I'm trying to find the DTFT of $$(-1)^n \cdot \frac{\sin (\pi n/2)}{\pi n}$$ I know the DTFT of $\frac{\sin \pi n/2}{\pi n}$ = a box function of amplitude 1, cutoff $\pi/2$. And I know that multiplication in time domain = convolution in frequency, but I'm still not sure how to piece everything together.


texts: Filtering Overlapping Frequencies in sound file

texts:
I have been given a sound file of a plane passing over a rain forest filled with birds. I am supposed to filter out the sound of the plane as it flies over.  I've accomplished this with various types of filters in MATLAB, but I always run into one problem. I can either cut out all of the plane and lose some of the background rain forest noise, or leave (almost) all of the rainforest noise and leave a decent amount of airplane noise. This is the frequency spectrum of the original sound file:

(Sorry for the crappy image.) Note the different colors represent left and right channels. From the plot, the majority of the plane's frequency is in the 0-1000 Hz range. The plane is clearly heard however until I use a high pass filter to remove everything from around 1700 Hz to 0 like so:

I suspect that the plane is causing that spike around 1700 Hz but evidently the rain forest is still producing lower frequency noise than that.

What would be the best method to remove that spike while keeping all that extra rain forest frequency?

EDIT:

Here is the plot of the time domain. Note how the plane gradually gets louder than fades away.

The sampling frequency is 22050 samples/sec.

And the spectrogram.


texts: Why do I have to duplicate the labels used for the classification process?

texts:
I'm using the BCI Competition dataset IV 2a. I'm working on Left & right hand movement.  I extracted data of channels (Fz,Cz,Pz).
I filtered data [8-30], then normalized it.
I calculated the common spatial patterns of the training dataset using the CSP.m of the Biosig toolbox
I want to use the "train_sc" to classify my data.
My question is why do I need to duplicate the labels of each event?
ex:
for each right hand event (750 samples) I must have 750 label = 1
for each left    hand event (750 samples) I must have 750 label = 2

labels = repmat(classlabel’,n,1)’;   %where n is the number of samples/event
CC = train_sc(features,labels,'LDA');

Does that mean that classification is done on samples rather than events?
Is each sample classified independently?!!!

The results I got are very disappointing. I've tried many classifiers but the best accuracy I got was 0.5
I need to know what I'm doing wrong.

I've read too much on CSP & classification of the eeg signals, but still I feel that I'm missing many things :(
Any help please.

Best Regards,
Sahar Selim


texts: Discrete Fourier transform - finding the fundamental quickly?

texts:
First, I appologize because I'm a software developper and it's been a very long time I didn't dive into pure mathematics, so my question might seem dumb. I hope not.

The context is pitch recognition in music.

If you take a musical note, and apply a Fourier transform to it, you will have and infinite sum of amplitudes for given frequencies. For example, if I play a note whose fundamental is $F$, on any instrument, after Fourier transform, I will have harmonics at at $F, 2F, 3F,\ldots,nF$. Every frequency will have a given amplitude which defines the timbre of the instrument (piano, voice, trumpet, ... all follow this loaw, but you'll have different amplitudes for every harmonic)

Now what I'd like to do is from a given audio signal, find $F$. Just that.
It's more complicated than it seems because you will always have background noise and so on ... Further more, $F$ isn't necessarilly the frequency with the highest amplitude !

So my idea for finding $F$ is to apply a DFT (well actually an FFT for speed) and find a frenquency $F$, so that $F + 2F +3F + \ldots + nF$ is maximal in the FFT output.

Do you think that's possible at all ?
Do you think that's possible in a very short time (let's say < 5 milliseconds) ?


texts: How to find location of minimum dimension of contour using opencv?

texts: suppose this is our input image and now we have to detect the point of minimum dimension, which I show as the red dot in an output image


texts:
I wonder do I understand beamforming right.
Question 1
Does same data inserted to all antenna at beamforming? Like figure below? The only different is through different phase shift? If it is right how about Massive MIMO which use 64,128 antennas? The all 64, 128 antennas get same data?

Question 2
Does concept of layer that I understand is right? Like figure below? I mean the two resource blocks transmitted at the same time and same frequency at same antenna?


texts: use of complex number in digital signal processing

texts:
how  complex number used in digital signal processing.
detailed information about on how this method will be useful in radar technology.


texts: How to avoid leakage in the following FFT

texts:
I have written a code for FFT for a signal which time domain plot and FFT is shown bellow.
I am not sure about what FFT length to do in order to avoid leakage of energy?
I know the number periods has to be integer, but i cant see how to implement this principle in here practically.
Also I get an amplitude which is very far from 5 5 10 as written in the time domain signal.
Where did i go wrong?

    clc
    clear all
    Fs=200e3;
    Ts=1/Fs;
    dt=0:Ts:5e-3-Ts;

    f1=1e3;
    f2=20e3;
    f3=30e3;

    y=5*sin(2*pi*f1*dt)+5*sin(2*pi*f2*dt)+5*sin(2*pi*f1*dt)+10*sin(2*pi*f3*dt);
    %plot(dt,y)
    nfft=length(y);
    nfft2=2.^nextpow2(nfft);
    fy=fft(y,nfft2);
    f_half=fy(1:(nfft2/2));
    xfft=Fs.*(0:nfft2/2-1)/nfft2;
    plot(xfft,abs(f_half));

fixed final:
clc
clear all
Fs=200e3;
Ts=1/Fs;
dt=0:Ts:5e-3-Ts;

f1=1e3;
f2=20e3;
f3=30e3;

y=5*sin(2*pi*f1*dt)+5*sin(2*pi*f2*dt)+10*sin(2*pi*f3*dt);
%plot(dt,y)
nfft=length(y);
nfft2=1000;
fy=fft(y,nfft2);
f_half=2*fy(1:(nfft2/2));
xfft=Fs.*(0:(nfft2/2)-1)/nfft2;
plot(xfft,abs(f_half)/(nfft));

Updated FFT conv FIR
clc
clear all
Fs=200e3;
Ts=1/Fs;
dt=0:Ts:5e-3-Ts;

f1=1e3;
f2=20e3;
f3=30e3;

y=5*sin(2*pi*f1*dt)+5*sin(2*pi*f2*dt)+10*sin(2*pi*f3*dt);
%plot(dt,y)
nfft=length(y);
nfft2=1000;
fy=fft(y,nfft2);
f_half=2*fy(1:(nfft2/2));
xfft=Fs.*(0:(nfft2/2)-1)/nfft2;
%plot(xfft,abs(f_half)/(nfft));

cut_off=25e3/(Fs/2);
order=64;
h=fir1(order,cut_off);
con=conv(y,h)
%plot(con)
fz=fft(con,nfft2);
f_half_z=2*fz(1:(nfft2/2));
plot(xfft,abs(f_half_z)/(nfft));

Updated with filter response plot:
clc
clear all
Fs=200e3;
Ts=1/Fs;
dt=0:Ts:5e-3-Ts;

f1=1e3;
f2=20e3;
f3=30e3;

y=5*sin(2*pi*f1*dt)+5*sin(2*pi*f2*dt)+10*sin(2*pi*f3*dt);
%plot(dt,y)
nfft=length(y);
nfft2=1000;
fy=fft(y,nfft2);
f_half=2*fy(1:(nfft2/2));
xfft=Fs.*(0:(nfft2/2)-1)/nfft2;
%plot(xfft,abs(f_half)/(nfft));

cut_off=25e3/(Fs/2);
order=32;
h=fir1(order,cut_off);
fh=fft(h,nfft2);
plot(abs(fh(1:nfft/2)));
con=conv(y,h)
%plot(con)
fz=fft(con,nfft2);
f_half_z=2*fz(1:(nfft2/2));
%plot(xfft,abs(f_half_z)/(nfft));
```


texts: How to take Processing/Scalloping losses in account in SNR/SFDR/THD/SINAD calculations?

texts:
I'm in the process of developing DSP library for our systems.
I have developed "Windowing" functions (Hanning, Hamming, Blackman-Harris, ...).
Each of these windowing function has his own "Processing Loss" factor.
For instance, in case of Hamming window, I have found that such Processing Loss factor is 1.34.
How should I take this factor into account when computing SNR, SINAD, SFDR, THD ?


texts: Can a linear reconstruction in compressive sensing perform well?

texts:
I am trying to implement compressive sensing for grayscale 2D images, then reconstructing them using a multi-layer perceptron(MLP). It seems to perform well no matter how many layers I add or remove, even without activation functions which is basically a linear combination of the inputs. there must me something wrong with my implementation. this is what I'm doing:
Encoding

Load image (PIL)
normalize image by dividing pixel values by 255
partition image into $B \times B$ blocks(patches) to get a matrix of size $P \times B \times B$ where $P$ is the number of blocks.
Flatten the blocks to form a matrix of size $P \times B^2$ let's call it $X$
Take M measurements from each block by getting the dot product of $\Phi$, an $M \times B^2$ Gaussian matrix, and $X^T$, the matrix containing the flattened, normalized blocks. The result is the measurement matrix $Y$ of size $P \times M$

Decoding
From what i understand, this is a sparse approximation problem. meaning my images(or blocks) are supposed to be sparse in some basis $\Psi$ in order to get an acceptable construction. I assumed the basis to be 2D-DCT. Therefore for each block, i try to reconstruct $S_i$ which is a $k$-sparse vector representing the flattened block in frequency domain. So my optimization problem becomes as follows:
$$\min|S_i|_1   s.t.   |\Phi\Psi^{-1} S_i-Y_i|_2 < \mu$$
$\mu$ being a small quantity. Now this applies to per-image optimization methods. What i am trying to do is an MLP which learns a universal function $f(Y_i)$ for reconstruction. So after training, my decoding process is as follows:  $$\hat{S_i} = f(Y_i)$$ following the acquisition of $\hat{S_i}$ i apply $\Psi^{-1}$ to it in order to get $\hat{X_i}$. finally, $\hat{X_i} * 255$ gives me the flattened block pixel values. reshaping each $\hat{X_i}$ into a $B \times B$ matrix then arranging them gives me the full image $X$.
Training
My dataset consists of $50,000$, $256 \times 256$ images. I just use the first batch (1000 images) since adding more data does not affect my model at all somehow. I do the aforementioned encoding process to every image then concatenate and shuffle all the blocks (i.e instead of images, i now have a matrix of flattened shuffled blocks from multiple images). My model inputs are the measurements taken from these images, of size $1000*P \times M$. To get my training labels, i apply 2D-DCT to each unflattened block before taking the measurements (right after encoding step 3 above), then flatten the sparse labels $1000*P \times B^2$, let's call the matrix $S$. Loss is calculated by taking $MSE(\hat{S}, S)$. Number of layers/nodes-per-layer do not seem to affect the model accuracy at all.
My questions

A Linear(no activations) MLP with a single hidden layer performs the same or sometimes better than a Multi-layer model with ReLU activations. What might be the cause of this?
In the literature, to guarantee construction $X$ has to be sparse in some domain. But i tried removing all transformations from my model and just worked in spatial domain and the model performs exactly the same. shouldn't it struggle to reconstruct the signal?

Output example
A model with 1 hidden layer having $B^2$ nodes takes the input of shape $M$. then the intermediate output is given to the output layer having $B^2$ nodes. in short, $M \rightarrow B^2 \rightarrow B^2$, no activations. training $MSE$ is $0.002$. $10$ of the blocks are reserved for validation, validation loss is also $0.002$. $B = 8$, M = $ 0.25 * B^2$. Now trying the trained model on an unseen images like "lena.bmp" $(512 \times 512)$ gives the following results compared to OMP:

N.B. The same $\Phi$ used in training is used in testing.
Any help is appreciated.


texts: Obtain the complex envelope of a two-tone signal

texts:
I know it is possible to represent a single-tone signal by its complex envelope:
\begin{align}
x(t) &= A(t)\cos(2\pi f_c t + \Phi(t))\\
x(t) &= \operatorname{Re}\left\{A(t)\mathrm{e}^{j \Phi(t)}\,\mathrm{e}^{j 2\pi f_c t}\right\}\\
x(t) &= \operatorname{Re}\left\{\tilde{x}(t)\,\mathrm{e}^{2\pi f_c t}\right\}
\end{align}

I think it should also be possible to represent a two-tone signal by its complex envelope:

$$x_2(t) = A_1(t)\cos\left(2\pi f_1 t + \Phi_1(t)\right) + A_2(t)\cos\left(2\pi f_2 t + \Phi_2(t)\right)$$

Can somebody show me the necessary steps to obtain the complex envelope of $x_2$ ?


texts: Why zero padding the 2-d DFT interpolates images in spatial domain?

texts:
I was applying different image interpolation techniques and I came know to about interpolation in frequency domain. In this technique we first take 2d DFT of an image, pad it with zeros and take the inverse DFT. The magnitude of inverse DFT gives interpolated image. It is working quite well!

Now my question is that how can it interpolate image in spatial domain when we are just padding zeros in frequency domain? How it constructs data points? And why only use zero padding?

Now I know about Shannon interpolation theorem (sinc interpolation) in which we can ideally reconstruct continuous signal from its discrete counterpart (if the sampling frequency is greater than Nyquist frequency). We can construct unknown values by sinc interpolation, but I can not understand how inverse DFT is constructing values for interpolated signal?


texts:
This is probably a very simple case but I am not sure where I am going wrong. In order to estimate the effect of path loss in Matlab, the function fspl(R,lambda) gives a pathloss value for a distance Rin $dB$.

My question is that, for using this function should one directly multiply the path loss value in $dB$ or convert it to linear scale and then multiply. If I use the $dB$ value, the signal power increases and if I use the linear scale value, there is a lot of difference between the transmitted and received power. Here's a snippet of my code:

SNR =2;
Fs=12e6;
Fc = 2.4e9;
d = 5; lambda = physconst('LightSpeed')/Fc;
data=randi([0 1],3000,1);
hMod = comm.OQPSKModulator('BitInput',true);
modulated_sgl = step(hMod, data);
pathloss_dB = fspl(d,lambda);
tx_sgl = modulated_sgl * pathloss_dB ;

If I try to plot the power spectral density of tx_sgl and modulated_sgl using the following:

hss = dsp.SpectrumAnalyzer('SampleRate', Fs);
step(hss,[tx_sgl modulated_sgl]);
release(hss);

I observe however that the power of tx_sgl is higher than the modulated_sgl. Can someone guide me where I may be going wrong.


texts: Finding $E_b/N_0$ when coding rate is $R_c$

texts:
I would like to calculate $E_b/N_0$  when the code rate is $R_c$ for BPSK signal. What I got this

$$
E_b/N_0= \frac{1}{2R_mR_c \sigma^2}
$$
and that is the case also according to this page.

$R_m=1$ for BPSK and $\sigma$ is standard deviation of the noise. Therefore I can say new $E_b/N_0$ is equal to old $E_b/N_0$ multiplied by $1/R_c$, therefore in $\log$ scale
$$
{E_bN_0}_{[\rm db]}={E_bN_0}_{[\rm db]}(\textrm{old})-10\log10(R_c)
$$
however this is different from what  MATLAB help says since according to them

CodedEbNo = UncodedEbNo + 10*log10(codeRate);

Can anyone explain where might be the problem?
And why my equation is not the same as Matlab help?


texts: How do I compare chords played on a guitar with a known sample?

texts:
I want to create a project that will recognized the chords that i played in guitar.
I recorded a raw sound to be compare to the one that the user will input.

Example

I have a recorded sound chordA then I will ask the user to play a sound in the guitar and check if the signal is similar to my sound chordA..

Getting the audio signal is easy but comparing them is so hard. What is the method that will suit my application?

This is my code:

recObj = audiorecorder;
disp('Start speaking.')
recordblocking(recObj, 4);
disp('End of Recording.');
myRecording = getaudiodata(recObj);
wavwrite(myRecording,'chordread.wav');
[x, Fs] = wavread('chordread.wav');
[x1, Fs1] = wavread('chordA.wav');
Y = fft(x);
Z = fft(x1);

Now that I have Y and Z how can I compare those two?


texts: Get spectrogram matrix without specgram in MatLab

texts: I have a signal in Malab, and I need a matrix, probably transposed specgram matrix without using specgram tool. Is it possible?


texts: Isolate g(t) starting from the definition of power spectrum

texts:
Starting from the following definition of power spectrum of a signal $g(t)$ in time:
$$S(f)=\int_{-\infty}^{+\infty}\rho(\tau)e^{-i2\pi f \tau} d\tau$$
where
$$\rho(\tau)=\int_{-\infty}^{+\infty}g(t)^{.}g(t+\tau)dt$$
Is it possible to isolate $g(t)$?


texts: Computer exercises and solutions in signal processing

texts:
I am very interested in this question and can literally not find anywhere a simple answer, so apologize in advance if it does not fit well this community.

I am looking for either an online resource or a book with computer exercises and solutions in signal processing (can be really anything in this space). I am interested in practicing independently, and while there are countless sources for practice computer exercises, I could not find any source that has solutions as well. Don't mind buying a $100 text book if needed. Appreciate any help in this direction.


texts: Spectrum analyser and resolution bandwidth

texts:
For the preliminary part of a lab I had to do, I was asked the following question:

I get how to do 1.1, but with 1.2, I really don't get what's going on; I really would've thought that the resolution bandwidth would be equal to the bandwidth of the low-pass filter.  Furthermore, I'm pretty sure a spectrum analyser has a variable band-pass filter in it, so I really don't get what's going on here.

Any help would be much appreciated.  For the record, this isn't part of an assignment or anything like that, it's not worth anything, it's just something that's bugging me.

Thanks.


texts:
I can only write .wav files to the current directory

fname = 'bassswoon.wav'
wav = wave.open(fname,'w')
wav.writeframes(struct.pack('h', bassswoon))
wav.close()

Adding a path to fname, i.e. fname = '/path/to/bassswoon.wav' does not work - the file is written neither in the intended directory ('/path/to/') nor in the working directory.

I don't see anything of use in the documentation - is it simply not possible with the wave object? Surely, there must be a way to hack it...


texts:
So I had this question (f) at my last exam that I failed and I'm now studying for my re-exam, the problem is that I still don't know what the answer is supposed to be in this question (f):

Do you have any suggestions how to solve it?


texts:
I am following the process that is described in this question: Transfer function of second order notch filter , I want to create a notch filter with the band suppressed equal to $f_c = 4000$ Hz, so using $\omega_n= f_c / f_s $, ($f_s = 48000$), I obtained the $\omega_n = \frac{\pi}{6}$, then using the exact same formula, with $a =0.8$. The pole-zero graph I get has the zero in $1$ and a pole in $0.8$, is this correct?? I am getting the half of the filter since the filter is centered in $0$ and not in $4000$ Hz.

As far as I know it must be centered in the wn I have (based on $f =4000$ Hz), but I am not sure why it is centered in $0$, or how to center it in the desired frequency.
I get a pole zero graph like this one, with zeros in $1$ and poles in $0.8$.

Please advice


texts:
I'm just learning about Fourier Transforms and as an input I'm using a WAV file (Matlab) and taking one channel of it and performing a quick Fourier transform. I've zero padded the input, as I've been told that it gives a better frequency resolution. I've done this and I see an identical looking graph apart from that the magnitude has changed. Could somebody please explain to me why the magnitude has decreased and is their any linearity to how it changes?

I assume that the frequencies will be similar to those before the zero padding but just be a better fit?


texts:
I try to describe a simple signal in the frequency domain with the use of FFT function:

from scipy.fftpack import fftfreq, rfft

def spectrum(sig, t):
    """
    :param sig: signal.
    :param t: time scale.
    :return:
    """
    f = fftfreq(sig.size, d=t[1]-t[0])
    y = rfft(sig)
    return f, y

Signal:

...and it's FFT:

Why I'm getting negative values of the FFT? Do I even use the FFT function well? I want to check "how many" frequencies are in the shown above signal. Thank you in advance.


texts: How would I calculate the state-variables of digital biquad filter (direct form II transposed) at steady-state?

texts:
I'm attempting to calculate the steady-state state variables for a digital biquad filter direct form II (transposed). Illustration

For example, let assume the filter is fed a constant input of magnitude C. What would be the best method to calculate the internal state variables as the filter approaches steady state?

When I run a simulation, it appears that S1 and S2 are negatives of one another (S1 = -S2). Is there a way to calculate the exact values?

My best guess is that I need to set-up a system of equations and then find the steady-state output (given constant input C). However, my calculations seem to fall apart when I actually try to doing them.

Edit: So after working this problem a little longer, I was able to determine the following:

$$
S_{1}[n] = X[n]\cdot (Gain_{DC} - b_{0})
$$

The DC gain is then found by taking the Z transform of difference equations and setting Z = 1.


texts:
I'm still new in OFDM and reading about it, but I have a question.

In OFDM system, when having $N$ data symbol assigned to sub-carriers, $N$ should be orthogonal with each others, which means $NN=0$. That's ok.  What's about that matrices $N$ itself?

Is there any property in OFDM saying that $N N^H$ should equal $I$ (where $I$ is the identity matrix)?


texts: I'm currently working on Person Identification, and using MFCC as the feature to my classifier. I have a database full of tagged(with labels) voices, and I use a CSV file to access the respective voice files. I generate the coefficients by dividing each voice file into equal frames of 1024 bytes & then, I find the coefficients of each frame. My problem is that all the training data isn't of the same size. Thus, the total number of coefficients vary from file to file. This dramatically affects the performance of my system. Is there any way around it? I've tried treating the whole file as a single frame, but the efficiency hasn't been very good. Any help would be appreciated.


texts: How to handle a logarithmic term in Kalman filter?

texts:
I am trying to implement a Kalman filter for an echo pulse detection application as similar to this paper. (an open source version is here (pg 16))
The measurement variable is $h(x,t)=A_0 (\dfrac{t-\tau}{T})^\alpha  \cdot \text{exp}{(\dfrac{t-\tau}{T})}$ and the state vector is $x=[A_0 \; \text{ } \alpha \; \text{ } T \text{ } \; \tau \text{ }]^T$.
In order to compute the measurement sensitivity matrix $H_k=\frac{\partial h}{\partial x}\Bigr|_{\substack{x=\hat{x}_{k^-}}}$, we need to compute the partial derivative with respect to $\alpha$. I obtain it as $\frac{\partial h}{\partial \alpha}=h \cdot log(\dfrac{t-\tau}{T})$.

$\tau$ would correspond to the start of the envelope of the pulse so it will be greater than first couple of time steps practically speaking. When I iterate the Kalman filter initialized with a real $\tau$ for example, then this partial derivative term becomes undefined for the start values of iterations because at $t=0, 1  \; time step, 2 \; time step \cdots$etc and so on are smaller than $\tau$ and the logarithm is undefined for negative values.


texts: Finding the Width of a Point Spread Function (PSF)

texts:
I have a binary image that has been convolved with the following PSF:

$$\mathrm{PSF}(x,y)=\cosh^{-2}\left(\frac{\sqrt{x^2+y^2}}{w}\right)$$ where $w$ is unknown.

The image binary image generally consists of contiguous black and white features that are larger than $w$. What is the most efficient way to estimate $w$? I typically don't have a problem to look at the image and provide an initial guess for $w$, and my goal is to get it accurately.


texts: Radar type application for narrow band only?

texts:
Not sure if this is the proper forum to ask this question, but still I am trying my luck.
I am looking for an application (in any area, medical, defence, or any other area) if exist, such that we want to dectect something by sending a waveform and get a retun on a linear array and determine the exact location of the target from the return signal. The catch is due to the situation or the physical limitation we cannot use wideband waveform. (We know that linear FM or wavefrom with higher bandwidth will help in spatial resolution, but we cannot use it due to some constraint.) In addition, the sampling rate cannot be increased so that we cannot form correlation matrix for better spatial resolution. This may also mean that target changes fast and we must rely only on time snapshot.
It appears in such case we stuck with Fourier beamforming.
The question is: Is there any such application where the above constaint really exists.
Any ideas will be highly appreciated.


texts: How Exactly Does MATLAB Zero Pad Signal?

texts:
I am doing some massive number crunching in MATLAB which involves millions of PSD estimations. Each data segment has length 41. So I have been using the multi-taper method with nfft=41 (details below in the code). My last simulation took more than eight hours. The code is optimized and vectorized as I can make it so I thought maybe if I change nfft=64 that might speed up. The profiler confirms that pmtm is the slowest part by far.

So, looking at the MATLAB help pages, I see that in the documentation it says nothing more than if nfft is larger than the length of the data segment then it is zeropadded. My question is what exactly is MATLAB doing? How exactly is the signal padded?

I can't really find anything googling and I ran my own experiments and I am not getting what I'm supposed to. Here I present my code and the result.

close all
clear all

dseg = [ 395.9640
  392.3630
  388.8230
  385.2820
  381.6130
  378.2600
  374.7470
  371.1860
  367.9920
  364.7830
  361.6630
  358.5660
  355.5620
  352.5820
  349.5530
  346.8110
  343.8130
  340.9010
  338.1460
  335.5110
  332.6300
  330.0470
  327.4020
  324.0930
  321.4680
  318.8760
  316.4470
  314.0850
  311.5310
  309.0600
  306.7570
  304.6690
  302.5820
  300.3890
  298.2480
  296.3650
  294.1660
  291.7610
  290.3300
  287.9700
  285.9850];

% Original signal with length=41
[ffty1,fftf1] = pmtm(dseg,4,41,1/30);
semilogy(fftf1,ffty1,'-o')
hold on

% Matlab "zeropads" to length=64
[ffty2,fftf2] = pmtm(dseg,4,64,1/30);
semilogy(fftf2,ffty2,'-or')

% Original zeropaded on both sides
dsegpad = [zeros(11,1); dseg; zeros(11,1); 0];
[ffty3,fftf3] = pmtm(dsegpad,4,64,1/30);
semilogy(fftf3,ffty3,'-ok')

% Original zeropadded at the end
dsegpad = [dseg; zeros(23,1)];
[ffty4,fftf4] = pmtm(dsegpad,4,64,1/30);
semilogy(fftf4,ffty4,'-og')

% Original signal just repeated
dsegpad = [dseg; dseg(1:23)];
[ffty5,fftf5] = pmtm(dsegpad,4,64,1/30);
semilogy(fftf5,ffty5,'-oc')

I compare the original PSD with MATLAB's padded PSD and they are close. Then I zeropad myself instead of letting MATLAB do it. I zeropadded on both sides of the signal. I tried zeropadding all on one side. I even tried periodically repeating the signal to pad the length. But in all three of the cases my padded PSDs are way larger than the original and MATLABs PSDs. Here is what I see.

I tried looking at pmtm code but got lost pretty fast. Anybody knows what MATLAB is doing and why does it obfuscate it so?


texts: Ideal low pass filter output at given sampling frequency

texts:
Consider the signal $\cos(30t)$ sampled at $w_s=40 rad/s$ using a unit impulse train. The sampled signal is filtered with an ideal low pass filter with unity gain and cutoff frequency $w_c = 40rad/s$. Find the resulting output.

The solution is given as $\frac{20}{\pi}(cos(30t)+cos(10t))$

Attempt:

Taking the Fourier transform of $cos(30t)$ gives $\frac{\pi}{2} (\delta(\omega-30)+\delta(\omega+30))$.

Using the unit impulse train:
$X(\omega)_{T_s} = \frac{1}{T_s}\sum X(\omega-k\omega_s)$ but this gives $\frac{1}{T_s}\sum \frac{\pi}{2} (\delta(\omega-70)+\delta(\omega+10))$.

If this is an acceptable approach, it's unclear where the cosine arguments and factor of $\frac{20}{\pi}$ is coming from. Help is appreciated.


texts: How to prove $\int_{\Omega} \sum_{i=1}^{N} f_i(x)dx$ is equivilant with $\sum_{i=1}^{N} \int_{\Omega} f_i(x)u_i(x)dx$

texts:
I have a 2D image defined on a region $\Omega$. Assume that the region can be separated into $N$ sub-regions $\Omega_i$ such that $$\forall i,j=1... N:\Omega_i \cap\Omega_j=\emptyset$$ and $$\bigcup_{n=1}^{N}\Omega_n =\Omega$$

Let $u_i(x)=1$ if $x \in \Omega_i$, otherwise $u_i(x)=0$. Then for a familiy $\{f_n:n\in 1..N\}$ of functions on $\Omega$ we have $$x \in \Omega_k \implies \sum_{i=1}^{N}f_i(x)u_i(x)=f_k(x)$$

I have a energy function that is defined as following
$$E_1=\int_{\Omega} \sum_{i=1}^{N} f_i(x)dx$$

Does the below equation $E_2$ equate with above energy function $E_1$? How to prove it?
$$E_2=\sum_{i=1}^{N} \int_{\Omega}  f_i(x)u_i(x)dx$$

That is my prove


texts: Usage of Digital Cosine Transform in lossy compression of images

texts:
I was reading about Digital Cosine Transform (DCT Type-I) and it's use in JPEG.
I have the following doubt. My teacher offered an explanation which wasn't very convincing.
What property of Digital Cosine Transform makes it attractive for compression of images? Why Digital Sine Tranform is not used?
More specifically, what is the intuition behind the usage of DCTs?


texts: What is the significance of a spike like autocorrelation function

texts:
I have some questions and doubts to which the answers are difficult to extract from text books. Hence I have posted this question here. Help would be extremely beneficial. Thank you.

(1) The normalized ACF of white noise is similar to delta impulse function. What is the purpose and benefit of this spike nature in signal processing and statistics?

(2) If the length of the spike for white noise is higher in comaprison to the length of the spike for a random process then is there a significance and what soes this convey? Does the length of the spike have a special name and a role in applications to signal processing?

(3) A delta function is not a random process. Then why is the ACF of white noise, which is closer to the delta function preferred?


texts: Relationship between "received" signal power and "transmitted" energy per bit?

texts: In digital communications textbooks (e.g., Digital Communications by Bernard Sklar), the relationship between the "received" signal power $P_r$ and the "transmitted" energy per bit $E_b$ is given as follows: $P_r=E_b R_b$ where $R_b$ is the bit rate. What confuses me in this equation is that the channel attenuation is not considered in the equation. Shouldn't the equation be as  $P_r=\alpha^2E_b R_b$ instead where $\alpha$ is the channel attenuation?


texts: Carrier sense multiple access, MAC layer

texts:
A medium access control MAC mechanism in wireless local area network and Ethernet is called Carrier Sense Multiple Access.

There are two different types of this CSMA

1- CSMA/Collision Detection

2- CSMA/ Collision Avoidance

Does anyone know what the advantages and disadvantges of each are?


texts: Other techniques like overlap save/overlap add

texts: I know overlap save and overlap add are used for long data sequence filtering. Are there any other similar or better techniques like these?


texts:
I have a stateful black box with four real-valued inputs and one real-valued output.  My problem is to predict the output at each moment in time, given the sequence of inputs seen up to that point.  During a learning phase, I can vary the inputs however I want and observe the output.  There's a little noise, of course, and the black box doesn't appear to be completely deterministic.

Specifically, I'm modeling a hard drive, and I want to predict the access time of the latest request given all previous requests.  I want a more black-box approach, though, because of the complexity of explicit models, and because I want this to work for other similar devices such as SSDs.

A couple of people have suggested that signal processing might be appropriate to analyze the sequences of input and output values.

Are there any ideas from signal processing that could help me to predict the output, or to characterize the input?


texts:
I'm posting this question to both check my own understand of DFTs (I'm relatively new to the subject), and to ask a few questions at the bottom regarding the DFT and its output. If anything isn't allowed then just let me know.

If we have a continuous signal in the time domain, $x(t) = cos(2\pi f t)$, we need to convert it to a discrete time signal.

If we sample this signal at a rate $f_s$,  which is the no. of samples taken per second, our equation becomes:
$$x[n] = cos(2\pi f n T_s)$$
$$\Rightarrow x[n] = cos(\frac{2\pi f n}{f_s})$$
where n is the sample number.

If we wanted to find the spectrum of this sampled signal, we could use the DTFT which would give us a continuous and periodic output. But because we're working with computers, the output must be finite, so instead we use the DTF which will give us a sampled version of the DTFT output. The DFT is given by:
$$X_k = \sum_{n=0}^{N-1} x_n . exp^{-j \frac{2\pi kn}{N}}$$
where N is the total number of samples, n is the sample of the input signal x[n], and k is the number of the frequency bin.

The frequency bin resolution is $\Delta f = \frac{f_s}{N}$, where N is the total no. of samples of x[n].

If we were to sample an input signal $x(t)$ with maximum frequency $f_{max}$, then by the Nyquist Sampling Theorem we want the sampling frequency $f_s = 2f_{max}$ to sample up to $f_{max}$ without aliasing.

I was wondering if my understanding of the DFT above is correct, and does this mean that, for each DFT bin, the DFT performs N summations and N multiplications to give the amplitude of the frequency at that DFT bin? And lastly, if I were to plot the amplitude vs frequency of the output of the DFT, would I only consider half of the spectrum (up to $\frac{f_s}{2}$) and then just double the amplitudes of the frequency bins below that?


texts:
Context: I'm trying use matlab to apply a single-pole filter to a time-domain ramp waveform that is generated by a sequence of time-shifted "RC steps" that are added together.

The time domain voltage waveform is

$V(t) = \sum_{k=0}^{N-1} V_{step}(1-e^{-(t-k\Delta t_{step})/\tau})u(t-k\Delta t_{step})$

where $V_{step}$ and $\Delta t_{step}$ are constants

Using the fft of a one-sided exponential decay, the unit-step, and the time-shift property of the fft I get a frequency domain representation of:

$\mathfrak{F}(V(t)) = V_{step}\left(\sum_{k=0}^{N-1}e^{-j\omega k\Delta t_{step}}\right)\left(\frac{1}{j\omega} +\pi\delta(\omega) - \frac{1}{\frac{1}{\tau}+j\omega}\right)$

So now I want to generate this complex-valued fft manually in matlab, multiply by a filter response, and inverse fft. To start, I'm checking the frequency response of the input and its ifft to make sure it looks right:

N=32;
Vstep=25e-3;
tstep=10e-12;
tau=5e-12;
ts=0.1e-12;

Nfft = 2^nextpow2(max([N*tstep/ts N*tau*5/ts])); %Get enough points for the whole ramp with at least 5 tau's per exponenetial
w=(0:Nfft-1)*2*pi/Nff;t
f=w/2/pi/ts;

timeshifts=sum(exp(-1i*(0:N-1)'*w*tstep));
step=[pi 1./(1i*w(2:end))];
expdecay=1./(1/tau+1i*w);

Vf=Vstep*timeshifts.*(step-expdecay);
vt=ifft(V);

So the frequency domain Vf looks reasonable in amplitude. However, when I take the inverse and plot vt, it is definitely not correct.

Where am I going wrong? I suspect there's something I'm missing with the fact that this is a DFT not fourier transform. Also I know there's some subtlety to the FT/DFT of the heaviside, particularly at w=0, and I know that the fft at 0 will just be the average (times N) and I'm not sure my script accomplishes this.

I do know that I can just start with the time domain and then fft it, but I'm rather curious now.


texts: How to do training of the Ant Colony used for classification?

texts: I'm using the Ant colony for my database of EEG signal for the detection of epileptic seizure. I read this article and I'm confused. Is this algorithm for training set to genertae rules or just for classification?


texts:
I want to create a microphone array that can be connected to my laptop. Previous discussions in stack exchange and other forums suggested me to get a USB audio interface. I would want to do a lot of real time signal processing on the sound arriving at the microphones (To be precise, I want to calculate the time delay of arrival of sound at the different microphones). But from what I saw, the audio interfaces are a bit expensive. But external sound cards are not! Here:

http://www.ebay.com/bhp/external-sound-card-for-laptop

Can I use an external sound card to receive input from multiple microphone (create an array, that is).

I will require 4 mics, but right now I need to do it with only 2 mics to test the algorithms. Sampling rate would be around 40Khz.


texts: Array bounds check o(1) (Integral image?)

texts:
Not sure if this is the correct place to ask but I have a question about techniques for processing an array of pixels with o(1) complexity. That is to get the sum of a subset of an array but if this subset extends beyond the bounds of an array taking those as zero. This is with a 2D array.

In reality the array can go up to [2000][2000] but lets take a 2D array that is sized [4][4] for simplicity.

array[4][4] = {0,  1,  2,  3
               4,  5,  6,  7
               8,  9,  10, 11
               12, 13, 14, 15}

Now let's say I want to get the sum of the array between [2][2] and [3][3], i.e. adding 10 + 11 + 14 + 15. This can be made in a o(1) solution by using a summed area table.

However I am little confused about how to keep this o(1) complexity providing:

the array subset extends outside the actual array
the array indices outside the actual array are taken as zero.
The algorithm computation is still o(1)

So for example lets say for this pixel sum I am asked to get the sum of [2][2] to [4][4]. This has a theoretical array of, numbers outside actual are denoted f (for fake):

array_theory[5][5] = {0,  1,  2,  3,  0f
                      4,  5,  6,  7,  0f
                      8,  9,  10, 11, 0f
                      12, 13, 14, 15, 0f
                      0f, 0f, 0f, 0f, 0f}

or in an image:

I think I have to use some sort of filter or image processing technique that I can't seem to find online or don't know the keyword and am flirting with it in my search results.

Any help would be greatly appreciated!


texts: Where is the "time and frequency" in a Discrete Wavelet Transform?

texts:
The DWT is motivated by the ability to see signal in time and frequency at the same time, if I understood that correctly.

I am using the following Matlab commands for some signal x:

[A,D] = dwt(x,'db4');

Then I plot A vs some vector t, representing 'time'. I do know the algorithm behind A and D: low (high) pass filter, then downsample by 2.

My questions: Where is the time axis and where is the frequency axis if I all what I have are the coefficients A (approximation) and D (detail)?

Second question: If all of this is a simple filtering and downsampling, I don't see what is new in DWT...


texts:
Is the Hilbert envelope (analytic signal) feasible to be calculated in audio frames/blocks in real-time, when using FFT with windowing.

What I'm specifically concerned with is that how does it not create or not be able to create discontinuities (in the end points) between successive frames?


texts:
So I am taking the Fast Fourier Transform of the following function:

$$ x[n] = \displaystyle\sum\limits_{i=0}^{5} A_{i} \cos\left(\frac{\omega_{i}}{\omega_{s}} n + \phi_{i}\right) $$

Where the sampling rate $\omega_{s} = 2 \pi 8000$, the frequencies $\omega$, the phases $\phi$ and the amplitudes $A$ are:
\begin{align}
\omega &= 2 \pi \{500, 522, 610, 675, 746, 825 \}\\
\phi &= \{0,0,0,0,0,0\}\\
A &= \{1, 1, -1, -1, 1, -1\}
\end{align}

The length of the signal is $N = 2048$ samples. No window is applied before taking the FFT.

The resulting FFT spectrum should have only a real component since the imaginary component should be zero. However, when taking the FFT and plotting the spectrum I get the following:

Why is the imaginary part not zero? It seems to be a discontinuity. However, I am not sure why it exists. Why does the imaginary part of the spectrum look like this?


texts:
The project I am doing is to detect stones in an image of a jewel, so as to get the idea of the deduction to be made in the gross weight of the jewel.
I have done image-segmentation in the HSV color space following this paper. If given a raw input like IMAGE1 below, the method I follow is three-fold:

using "gold" pixels threshold values I detect the gold part in the image (IMAGE2 below)
then I detect the background/foreground in the image as the background is known to be white (IMAGE3 below)
intersection of the outputs of the first 2 steps -- not GOLD && FOREGROUND => STONE (IMAGE4 below)

The steps are illustrated in this image:

Some sample images and their corresponding results can be found:

image-1 -- result-1
image-2 -- result-2
image-3 -- result-3
image-4 -- result-4

Apart from this, I have tried using the flood-fill algorithm to flood the gold / stone pixels based on some pre-known seed points. But due to the non-uniform illumination and the fineness of the image, flood-fill didn't produce good results. Also the lo_diff and up_diff values (in the OpenCV implementation) are not universally fixed and need to be changed for each query point.

I also tried using template matching in it's naive form but the complications of variance of the occurrence of the template in the query image led to that method being not applicable everytime. And if the stone type is irregular, then template matching gives erratic results.

Segmentation in the HSV color space yielded better and more generally applicable results than the previous two methods. But as visible in the results of the 4 sample images previously mentioned, that too has it's own limitations, like

In the second image, due to poor illumination all of the background is (wrongly) detected as being foreground.
In the first image, due to (very) slight illumination problem in the left-top corner, the detected area is blown up.
Similar noisy issues are visible in the last 2 results.

Is there a better way of solving the problem? Can someone suggest some improvements to the current approach?
Suggestions at an altogether new approach are also welcome. Is it possible to implement a machine learning approach to the problem apart from image-processing? Can someone give any specific pointers in that direction?


texts: Modelling Unwanted Signal in a LMS Adaptive Filter

texts:
I'm having some confusion learning about the LMS Adaptive Filter. I understand that the whole model of adaptive filters relies on the fact that we give it a reference signal to which it keeps comparing the input * filter with and the filter coefficients keep changing until the error between input and reference is zero.

Let's say I have a telephone conversation with a time-varying sinusoid added on top that changes frequencies every so often. How exactly do I give it a reference signal? I mean, if I had a reference signal of what I wanted the input to become, I would have just used it? And similarily if I had a reference signal for my unwanted sinusoid, I could have just subtracted it from my original sound.

What am I missing here?


texts: what is the main differences between the methods for designing FIR filters?

texts:
what is the main differences between the methods for designing FIR filters ?
I mean specifically these methods

fourier transform method
window method
frequency sampling method
optimal design method (remez algorithm)


texts: Baby Cry recognition robust against the distance

texts:
I am building a baby cry detection system using deep learning algorithm, Raspberry Pi and a microphone as the sound sensor.

When the baby starts crying then the system should detect the cry. It is working good when I placed the mic near to the baby and but when the mic is away I am not able to detect sound.

Please suggest me any techniques to detect sound from max 3-meter range.


texts: Detecting phase difference between a clean and a noisy signal

texts:
I have a clean square-wave signal and a noisy signal from an accelerometer,
perfectly synchronized to each-other but out of phase. I want to find the phase difference between them in real time, in software (I am building a dynamic balancer).
I plan to convert the clean signal to a (complex) phasor in real time and multiply that by the noisy signal, then integrate over a moving window. I should be able to measure the phase angle of the integral to get the answer.
Does this sound like the right approach? It would seem to cancel out the other frequencies automatically.
I have considered measuring the frequency of the clean signal and filtering the noisy signal using a band-pass filter around that frequency, but it seems unnecessary.


texts: DSP Demodulation advice sought - please help

texts:
my field of engineering is condition monitoring. I use vibration analysis to determine the current state of health of a machine and then recommend any course of action required to ensure it does not fail during a production period. See www.efftek-uk.com for further description.

The most important tool I have is ESP which stands for envelope signal processing in other words, demodulation. The instrument I use to collect my vibration readings does this in the analog realm. The steps involved in this are as follows. The signal from the accelerometer is first passed through a high pass filter, then it is rectified and finally passed through a low pass filter before acquisition. The values of the LP and HP filters are selected according to the speed of the machine being tested for instance 2.5 & 5 kHz filters.

Another point which may not be relevant is that in VA, we use spectras with multiples of 400 lines, the top 112, 224 etc. are discarded for anti-aliasing. For example, we would use a 2.5-5kHz ESP filter and take 1024 samples at 1280 samples per second giving us a 512 line spectra of 0-640Hz and a final 400 line spectra of 0-500Hz

I have developed an ADC board for the Raspberry Pi (available of eBay at £19.50 - shameless advertising!!!) which was originally for use with displacement sensors so I could plot shaft centreline positions during operation of a gearbox to detect bending forces. I have decided to experiment using the board with MEMS accelerometers and would like to try and emulate what my very expensive data collector does in the analog realm with DSP. This is where I would like some advice.

I have already written code to sample, window and do fft's so I can look at the unfiltered acceleration and velocity spectras. Now I would like to implement filtering especially demod filtering.

I did read a thesis paper a while ago from a Swedish student who did demodulation using matlab and his results were encouraging but a) I don't have matlab and b) I can't find his paper.

I don't want anyone to write my code for me (but I may be back if my code doesn't seem to be living up to expectations) I simply want some advice on the process involved in DSP demodulation - is it the same as I mentioned above - HP filter, turn all values positive (rectification) then LP filter followed by FFT? and if so, what filters do I use - as a newbie to DSP, I don't know the difference between IIR and FIR etc. - what should I use?

I am coding on the Pi in C++ and have installed the Aquila-DSP library though as of yet, I haven't made use of it, I wrote the window function myself as the only window recommended for VA is Hanning and that was quite straightforward.

Thanks for reading, any advice will be gratefully accepted and do go to eBay if you are interested in an 8 channel 12 bit ADC with trimmers so the input voltage can be calibrated for each channel up to 0-20V for the Raspberry Pi.

Best wishes, Steve.
s.ferry @ efftek-uk.com


texts: Allpass group delay formula and effect on Karplus-Strong

texts:
Some introductory questions:

What is the unit of the group delay in a digital filter (is it samples?)
What is $\Omega$ in the equations below (is it $\Omega=2 \pi\frac{f}{f_{sample}} ?) $

If I got any of these wrong, then there is no point reading the rest of this question. Just correct me.

What I want to achieve

I am using a Karplus-Strong algorithm to simulate a steel string. I noticed, that I need to detune the harmonics to higher frequencies, so it really sounds like a steel string (see here). To achive this, I added an allpass to the delay in the feedback loop.

$$H(z) = \frac{a + z^{-1}}{1 + az^{-1}}$$

However the allpass adds an additional delay and detunes my string (lower). So I tried compensating for that additional delay by subtracting the group delay of the allpass from the main delay.

The result was: it still detunes, but in the opposite direction.

The allpass delay:

Matt L.  gave me this formula for the group delay

$$\tau_g(\Omega) = \frac{1 - a^2}{1 + 2a \cos(\Omega) + a^2} $$

$$\Delta t = \frac{\tau_g}{f_{sample}}$$

But robert pointed out that it is the phase delay which determines the pitch and gave me this formula:

$$ \tau_\phi(\Omega) \ = \ \frac{1}{\Omega} \arctan\left(\frac{(1-a^2)\sin(\Omega)}{(1+a^2)\cos(\Omega) + 2a} \right) $$

But alas, when I subtract either of the two from the total delay in the feedback loop (i.e. the regular delay plus the allpass) it appears I am overcompensating. There is no problem with posive $a$ which leaves the fundamental undelayed, but for negative $a$ (e.g. $a=-0.93$) the string is about half a semitone too high.

Edit (1)

I had an error in the formulas which is fixed by now (thanks to Matt L.). Originally my string was tuned too low, after the fix it is too high.

Edit(2)

Meanwhile (thanks to Matt L. and robert) I have little doubt, that my calculations are basically correct. I ran my allpass through scilab and the group delay is reported as apx. 30 samples. With the formulas from Matt and Robert I get values around 31 samples. Hoever, to tune the string correctly (by ear) I must compensate for 24 samples only.

Interestingly Miguel Negrão who did some research on plucked steel strings wrote in his thesis: Generally a negative value of a would have been used, which would have delayed the lower frequencies more than the high frequencies (phase delay), and would require a re-tuning of the main delay in the string loop to keep the string at the same fundamental frequency. In the SupaStrings system case, after some listening tests, a positive value of
a was used instead (non-physical) which keeps the low partials more a less in the same frequency but detunes the high partials down to a lower frequency.

It seems like he faced similar difficulties and couldn't solve them, even though he is quite a brilliant guy.


texts: Getting the bitstream of an unknown RF signal

texts:
I am trying to get the bitstream of a signal sent over 433MHz. I have installed SDR# (I'm on Windows), connected my NESDR SMART, and I identified the signal here:

I can see there is only one peak with progressively decreasing peaks, so I assume this is AM. Then I set the bandwidth to the width of the central lobe. I finally made a recording of the "Audio" and opened the WAV into Audacity but I could not make sense of the sines (I was expecting a digital stream, since the modulation was specified) so here is the I/Q data instead.

This is an extract of the body of this waveform:

I have also installed DSD+ and the SDR# plugin as well as the virtual cable between the two, but after setting up SDR# the waveform does not budge so nothing comes up.

What sort of encoding is this? How can I read the bitstream? I am open to other things to try to get to the bottom of this.


texts: How to correctly pass filter states from one processing block to the next

texts:
I'm processing a 4 channel audio signal and due to memory requirements need to process the signal in an arbritrary number of blocks. To ensure consistent filtering between blocks I've been told I need to pass the filter state at the end of the current block, into the filtering instance for the next block.
I'm using scipy.signal.lfilter() and this is used with an array of FIR filter coeffs and a loop to filter the signal into the number of frequency bands required for analysis. So I end up with an array of, for examples, 20 filtered copies of the original signal and dimensons freq_bands x time_step x channel. The original method for this (without passing filterstates) was using list comprehension seen below
filt_audio = np.array([signal.lfilter(filt,[1.0], audio_np, axis=0) for filt in filters])
This however doesn't work when trying to also generate and store filterstates.
On the first block, I create initial filter states for each channel using
        if zi is not None:
            zi = zi
        else:
            zi = np.array([signal.lfilter_zi(filt,[1.0])*audio[0,channel] for channel in range(audio.shape[1])]).T

I then attempt to filter the original audio signal using my array of FIR coeffs and generate/use the filterstates using the below loop
for indx, filt in enumerate(filters):

            filt_audio_temp, zi_temp = signal.lfilter(filt,[1.0], audio, axis =1, zi=zi)
            filt_temp.append(filt_audio_temp)
            zi_temp_lst.append(zi_temp)

        filt_audio = np.array(filt_temp)
        zi = np.array(zi_temp_lst)

However attempting to then filter the signal using my filterbank and the initial state for each channel gives the below error.
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-36-f3a7027a0591> in <module>
      4     block = audio_np[i*blocksize:(i+1)*blocksize]
      5     print(i)
----> 6     azi1, elev1, psi1, zi, past_sample = extractSpatialFeatures(block, fs, filt_taps, first_band, channels, rho = 1.2041, Z_0 = 413.3, c= 343, zi=zi, past_sample = past_sample)
      7     azi_temp.append(azi1)
      8     elev_temp.append(elev1)

<ipython-input-35-e18635346eb2> in extractSpatialFeatures(audio, fs, filt_taps, first_band, channels, rho, Z_0, c, zi, past_sample)
     26         for indx, filt in enumerate(filters):
     27
---> 28             filt_audio_temp, zi_temp = signal.lfilter(filt,[1.0], audio, axis =1, zi=zi)
     29             filt_temp.append(filt_audio_temp)
     30             zi_temp_lst.append(zi_temp)

~\anaconda3\lib\site-packages\scipy\signal\signaltools.py in lfilter(b, a, x, axis, zi)
   1850                         raise ValueError('Unexpected shape for zi: expected '
   1851                                          '%s, found %s.' %
-> 1852                                          (expected_shape, zi.shape))
   1853                 zi = np.lib.stride_tricks.as_strided(zi, expected_shape,
   1854                                                      strides)

ValueError: Unexpected shape for zi: expected (144000, 127), found (127, 4).

This seems to be saying it expects a filterstate for each timestep of the signal (144000, 127) however what it is receiving is the filterstate for each channel (127,4). 127 is the number of filter coeef. I also see zi needs to be zi.T. My thinking would have been that I just need a filterstate for each channel for the first block and then for the second block I'd need to pass states for each frequency band and each channel. So I would need a zi matrix of freq_band  x channel == 20 x 4


texts: Inverse DFT: Is there a valid / intuitive interpretation of results for non-integral timestamps?

texts:
I implemented the plain DTF / inverse DFT algorithm in C++ in order to help me understand the method. As a sample input I considered
f(x) = sin( x * PI / 5 )
and collected samples for x = 0, ..., 9. So my input data was (0, 0.58778, 0.951057, ...).
Computing the DFT formula gives two non-zero coefficients: c_1 = (0 - 5i) and  c_9 = (0 + 5i).
Now, I compute the inverse DFT. For x=1, these coefficients multiplied with their basis functions give
(1/10) * exp( 2 * PI * i * 1/10 ) * (0 - 5i) = (0.293893, -0.404508i)
and
(1/10) * exp( 2 * PI * i * 9/10 ) * (0 + 5i) = (0.293893, 0.404508i)
respectively which is correct as the real parts sum up to 0.58778 (thus, reconstructing the input value) and the imaginary parts sum up to 0.  This works great for integral timestamps x=0, 1, 2, ...
My question: ... but how to interpret the results if I plug in x=0.5 (result 0.000000 - 0.951057i) or x=1.5 (result 0.000000 - 0.587785i) for example? I see that it is not DFT's promise to provide a good interpolation here. But is there an intuitive interpretation?  (Came across this because I wanted to plot the DFT/IDFT result as a time-continuous graph and don't know how to deal with the complex results between)
EDIT
I add my C++ code here:
void dft( const std::vector<double> & vInput, std::vector< std::complex<double> > & vCoeff )
{
    const std::complex<double> i( 0., 1. );

    const int nElem = vInput.size();

    vCoeff = std::vector< std::complex<double> >( nElem, std::complex<double>(0., 0.) );

    for( int k = 0; k < nElem; ++k )
    {
        for( int j = 0; j < nElem; ++j )
        {
            vCoeff[k] += exp( -2. * M_PI * i * std::complex<double>( j, 0. ) * std::complex<double>( (double)k / nElem, 0. ) ) * vInput[j];
        }
    }
}

std::complex<double> inverseDft( const std::vector< std::complex<double> > & vCoeffs, const double t )
{
    std::complex<double> result( 0., 0. );
    const std::complex<double> i( 0, 1 );
    const int nElem = vCoeffs.size();

    for( int j = 0; j < nElem; ++j )
    {
        result += vCoeffs[j] * exp( 2. * M_PI * i * std::complex<double>( j, 0. ) * std::complex<double>( t / nElem, 0. ) );
    }

    return (result / nElem);
}

This is how I define my input:
std::vector<double> vInput;

for( int i = 0; i < 10; ++i )
{
    vInput.push_back ( sin( i * M_PI / 5 ) );
}

And this is the output I did not expect for t=0.5, t=1.5, etc.:
inverseDft( vCoeffs, 0. );  // okay (0.000000, -0.000000)
inverseDft( vCoeffs, 0.5 ); // ?? (0.000000, -0.951057)
inverseDft( vCoeffs, 1. );  // okay (0.587785, -0.000000)
inverseDft( vCoeffs, 1.5 ); // ?? (0.000000, -0.587785)

EDIT2
Great to have these expert answers here!  I can contribute the 3D plot which supports my intuition and therefore might also be helpful for others. It shows f(t) computed by IDFT and how it travels over the original data points (crosses). In order to achieve this, it makes use of the whole space including the complex plane.


texts:
I have an arbitrary shape defined by a binary mask (gray = shape, black = background).

I would like to find a largest possible rectangle containing only gray pixels (such rectangle is pictured in yellow):

The shape is always "one piece" but it is not necessarily convex (not all point pairs on the shape's boundary can be connected by a straight line going through the shape).

Sometimes many of such "maximum rectangles" exist and then further constrains can be introduced, such as:

Taking the rectangle with its center nearest to shape's center of mass (or center of image)
Taking rectangle with aspect ratio nearest to a predefined ratio (i.e. 4:3)

My first thought about the algorithm is the following:

Compute distance transform of the shape and find its center of mass
Grow square area while it contains only shape's pixels
Grow the rectangle (originally a square) in width or height while it contains only shape's pixels.

However, I think such algorithm would be slow and would not lead to optimal solution.

Any suggestions?


texts:
I have a simple program that captures audio from an audio device. I have configured a nominal sample rate of 48000Hz and a buffer size of 1 millisecond. The audio sub system should execute my capture code with approximately 1000Hz.

My assumption is that the audio clock (LRCLK) is not running at exactly 48000Hz. Using this 1kHz-callback from the audio system I want to approximate the actual sampling frequency relativ to another clock. Currently I use the system clock of my workstation as a reference. I plan to use a network clock.

For my very first experiment I used a simple moving average with a period of 1000 but there is still very much jitter. Sometimes the average jumps by 2Hz which seems much to me.

What is the proper way to go here?

Many thanks!


texts: From pressure to SPL : microphone calibration using sound level meter calibrator and sound level meter

texts:
I have a microphone + signal conditioner with adjustable gain and filter. I capture amplified signal using oscilloscope. I'm using this signal to compute equivalent continuous sound level

With signal conditioner, I set needed gain and A-weighted filter.
From oscilloscope signal V I can get pressure p = V/ (gain x mic sensitivity)

We can suppose A-weighted filter filters out frequency out of audio range. So I can easily compute leq = 20 log10(rms(pa)/pref), pa is time domain pressure A filtered, pref is reference pressure = 20e-6 pa

I used a sound level calibrator to check if it works. I get a leq of ~94dB and 114dB when applying 94dB and 114dB level with calibrator.

I also have a sound level meter, I also checked it with sound level calibrator.
Now, in my room (some noise because of oscilloscope...) with A filter it measures ~46 dB(A). Now I put my microphone close to it, with microphone I get a leq of 55 dB. When changing sound level meter filter from A to C results are better : sound level meter measures ~56dB.
I changed microphone acquisition duration but leq results are same.

How can we explain such leq difference when sound level meter and
microphone measure dB(A)?
How can I calibrate my microphone with    sound level meter?


texts: Segment/Idenfication of signal from a spectrogram

texts:
I have the following Spectrogram of a signal:

Essentially, in this signal, I know there is one point of interest (between 60-80) and between the the particular frequencies.

I could, potentially smooth the signal and perform a simple peakfinding algorithm to find the peak and then just extract the bin number where the bin has the highest frequency.

I'm however looking for a different method where multiple calls, at different bins of the output, for example:

Here there are 4 calls or 4 significant components of this signal, that I need to extract and show in a plot (not a spectrogram). Ideally, in this spectrogram, the X-axis represents time, so, I would like to extract the signal so that it is possible to show where they are in time.

The idea behind this, is to be able to count the number of calls (significant parts) of the signal.

I have seen this being implemented through the idea of Principle component analysis, however, I do not want to go down this route, instead, would like to use a simple and more elegant solution.

Could someone please recommend an algorithm, or, some reading materials to where I might be able to implement such an algorithm to be able to count the number of calls?

Thanks


texts: For extracting the characteristics of the signal what minimum duration should be considered?

texts:
Being non signal processing student I have a limited understanding of concepts.

I have a continuous periodic bearing faulty signal of 10 seconds duration sampled at 12 kHz. So I have 120,000 discrete time samples in 10 sec signal. I have utilized some machine learning techniques (Convolution Neural Network) to classify faulty signals to the non faulty signals.Based on RPM, one shaft cycle need 400 time samples approximately. Some research papers have considered the signal by segmenting by 1024 window size (non overlapping) to create 117 segments from 10 sec signal which are utilized as input to the network. So here in each of the 117 segments have approximately 2.553 shaft cycle (1024/400) information. I am confused because but no one mentioned why they are segmenting for such duration

My question:
Is the information from one shaft cycle enough to provide all details ? Because number of discrete time samples in one shaft cycle are always below the sampling frequency. If I segment it for 10 shaft cycles (for example), is there any explanation for what difference will it make?

Normally to analyse what is inside the signal are there any minimum duration's standards defined ?


texts: matching the mouth shape of an audio file(lip sync)

texts:
I am working on something like cartoon, but the character's mouth shape is automatically detected and matched by the word they are saying.

my guess is to use speech recognize algorithm to analyze the wav file, and get the content(word) at the exactly moment, for the exactly word, prepare corresponding mouth shape picture.

but it look like a dumb way... is there any convenient way to work this out?

I saw many games support this, how did they achieve that?


texts: How to flatten this OCT image with respect to the marked line

texts:
I am working on a medical image processing project where I have to use retinal OCT images. I am using the dataset available here. This dataset has a lot of images and I want to flatten all of them. I have a flattening code which works correctly only for some images but fails as a generalised solution to all the provided images:

For instance, if the input image is:

The output after flattening is:

The code that I had used for performing the flattening process is:

function img_shift = img_flatten(img)
%img = rgb2gray(imread('retina.png'));

    % Find bottom line indices
    %find the location of the curve which is white/gray in color different
    %from the black
    threshold = 100;
    retina = NaN(1, size(img,2));
    for j = 1:size(img,2)
        retina(j) = find(img(:,j) >= threshold, 1, 'last');
    end

    % Fit polynomial
    x = 1:length(retina);
    p = polyfit(x, retina, 2);
    retina_fit = round(polyval(p, x));

    % Plot the fit of the retina
    %figure(1);
    %plot([retina', retina_fit']);

    % Shift the data in the image
    img_shift = zeros(size(img), 'uint8');

    for j = 1:length(retina_fit)
     img_shift(:,j) = circshift(img(:,j), -retina_fit(j));
    end%for j

    % Plot the image
    %figure(2);
   % subplot(2,1,1);
    %imshow(img,[]);
    %subplot(2,1,2);
    imshow(img_shift);
end

I tried changing the threshold but it didn't make any difference to the code.

I have the following questions regarding my code:

1) Is there any other better method to perform flattening of OCT images?

2) Can I tweak my code to have a generalized flattening code for all the images belonging to the dataset?


texts: Decimation and filtering in the frequency domain

texts:
I'm working on a Software Defined Radio project where I'd like to low-pass filter and decimate an analytical signal (IQ) sampled at 96ksps. Let's say the low-pass filter has a cutoff at 5kHz and I'd like to decimate by a factor 4 so that I have 24ksps out.

The idea is to perform the filtering using fast convolution using the overlap and save method as described in this article [pdf]:

http://www.3db-labs.com/01598092_MultibandFilterbank.pdf

I'm wondering if there are any pitfalls to my approach:

Performing an N length FFT.
Then doing an N length circular convolution (my multiplying with the FFT of my filter of length P).
Then performing an N/4 IFFT back to decimate by 4 using the N/4 center taps of the forward FFT. Since my filter is a low-pass with a cutoff at 5kHz there should be very little energy outside the N/4 center taps of the FFT, and the P - 1 samples I need to discard should also lie outside the IFFT (if my filter is not too long).

EDIT:

This specific application is on a Raspberry Pi 3. After having given this some more thought I've realized it's not as clever as I first thought. I've stared too much on 0 centered FFTs and briefly forgot this is not the case. I would have to "remove" zeros in the middle of my FFT to make it shorter and then perform the IFFT.

What I will do is to do the fast convolution with an N length FFT and N length IFFT, and decimate when I copy samples to the output buffer.


texts: Poles of the transfer function in Z Transform

texts:
Given that:

$H(z)$ has 4 poles maximum.
$H(z)$ has a pole at $z_1=a+bi$

Given that the impulse response $h[n]$ is:

Symmetric: $h[n] = h[-n]$
Real: $\forall$$n$ , $h[n]$$\in$$\mathbb{R}$

How we can conclude that the other poles are the inverse and the complex reflection of $a+bi$ ?


texts: Can you downsample unevenly with even distribution?

texts:
Could I downsample a series of length 10 to length 7, with even distribution?
i.e. Numbers get split across two (neighbouring) resultant numbers, as opposed to one result being a combination of 1 number, and another result being a combination of 2.
i.e. All numbers be a combination of 10/7 numbers.
For this example, with my current - naive - approach, some of the 7 "buckets" would average 1 number, and some other buckets would average 2 numbers: avg(1), avg(2, 3), avg(4), avg(5, 6), avg(7), avg(8, 9), avg(10)
Would the solution be something like: weighted_avg((1, 1), (2, 10/7 - 1)), weighted_avg((2, 1), (3, 10/7 - 1)), ...
where weighted_avg = ((value_x_1, weighting_1), (value_x_2, weighting_2)) -> weighted average


texts: Given the input and the output, how to determine the impulse response?

texts:
I would like to find the impulse response, $h[n]$, of an LTI system given the input

$$x[n] = [1,-3,2]$$

and the output

$$y[n] = [1,-1,-4,4]$$

I know that $y[t]=x[t]*h[t]$, but I am having hard time to figure out the right way to calculate the impulse response. I know very little about signal processing, so if you don't mind giving an easy explanation, then I appreciate it. Or, if it's possible to do an example, that's better.


texts: Given an FIR filter coefficients, how can I find coefficients of its equalizer filter?

texts: I have an N-tap direct form complex FIR filter with known coefficients [c0 c1 .. cn]. I want to find the coefficients of its equalizer filter which when cascaded with the original filter gives constant magnitude response for all frequencies. How can I do that? Thanks.


texts:
I'm a bit confused with some details about frequency domain or studying systems with Laplace and Fourier analysis.

Let's say we have a system with a transfer function $H(s)$ with no poles and zero's at origin. We know that the transfer function is a ratio between output and input. Also, an input of $X(s)=1$ represents a single impulse (Dirac) in time domain, which can be viewed as a signal containing all frequencies in frequency domain. Also, applying this input to the system will lead to an output $Y(s)$ which is equal to $H(s)$ which is like the system's unique characteristic (its transfer function) and also can be viewed as a response to all frequencies. This signal $H(s)$ can be transformed back to time domain.

Let's say we have a system whose impulse response is something like this:

When I studied Fourier Transform we used to study system's response and Bode plot by analysing variations in $\omega$. So, for example, making $\omega=0$ we could find the system's DC gain. Making $\omega$ low values we can observe the system response to low frequencies and the same for high frequencies with high values in $\omega$. It was very useful to study filters, for example.

But note that all of this is only using the transfer function. And the transfer function may have an appearance like this above. And this does not make sense to me. If we think of the time->frequency transform, we should think that in frequency domain we are finding all frequencies which compose our signal in time domain, correct? This explains the Dirac containing all frequencies, for example. So, when making s=0 we should find H(s)=0 since this signal above does not have a DC level.

When we think about the transfer function as a system response to $X(s)=1$ it seems ok. But if we think about the transfer function as the same representation of this signal above, it does not make any sense. Also, if we apply a step signal, our system may respond like this:

We saw that for this system (with no poles at origin), we can find something like $H(0)=K$ which represents the system's DC gain. Applying a step response will lead us to a signal $\frac{H(s)}{s}$. So now making now s=0 will lead us to $+\infty$. Again, if we think of this signal being the one above but represented in frequency domain, making s=0 should left us a DC finit level since the time domain representation show us that the signal contain a limited DC level.

I know the final and initial value theorems but still I want to understand these relations between these domains.

Note that if we want to study a filter response, we should apply some signal to the input and see what is the output. Mathematically we should convolute these signals in time domain or multiply in frequency domain and see what is the resulting signal. In time domain it all makes sense since we just convolute and see the system's response. In frequency domain it also makes sense just like the case where you use $X(s)=1$ which represents all frequencies and multiply by the system and see the resulting signal. But when thinking about the same signal being represented in both domains, I think I'm messing things up.

I just want a help to understand relations between 's' variable and time domain. Making s=0 shouldn't lead me to the DC component of the signal in time domain? Generically, each value of 's' shouldn't lead me to 'how much' of this frequency is present in this signal?

in short:

A signal like the first one does not have any DC component so it goes to zero when $t=\infty$. So why its representation in frequency domain have a DC component when we put s=0?
Also, the second signal does have a DC level when $t=\infty$. But in frequency domain the signal goes to $\infty$ when s=0 (which I thought it should represent DC level).
It happens with Heaviside function too. Its representation in frequency domain is $\frac{1}{s}$ altough when s=0 it goes to $\infty$. But it has a DC limited component. How do we explain all of that?


texts: Oversampled Binary Image Sensor vs Active Pixel Sensor(CMOS)

texts:
I am a beginner in the field of Image processing.I was designing a camera sensor when I came across these two type of sensors:

CMOS(Active Pixel Sensor)
Oversampled Binary Image Sensor.

Now, I can find many websites for comparing CMOS vs CCD but nowhere can I find a comparison with OBIS. Why is this? Is the OBIS too outdated to be used? I was thinking of designing a camera that captures images using OBIS and later converts them into RGB format (with a variety of processes I have not discussed here). As advantages of Binary Sensor I mentioned Low image size and super fast processing. Are these really going to be an advantage over CMOS sensor? Will my camera be a success or is this a really bad idea trying to use Binary Image sensor? Is using Binary Sensor a good idea or is it a too outdated idea?

Thanks in advance for helping this novice!!


texts: Radix-8 butterfly with Winograd and Cooley-Tukey algorithm

texts:
I saw the Winograd radix-8 kernel algorithm below, shown in the image. Comparing to the mathematical formula of Cooley-Tukey, there is a multiplication by $\cos$ and $\sin(\pi/8)$, which can't be easily realized by the combinations of $1$ and $\sqrt{1/2}$, which are the components used in Winograd.

As we can't realize something out of nothing, what is the mathematical trick to prove the Winograd and Cooley-Tukey are equivalent? Are there any such proof online I can read?

Winograd algorithm:

Cooley-Tukey algorithm: