texts: Nonnegative or positive band-limited interpolation

texts:
Given samples of an everywhere non-negative or positive-valued continuous-time signal band-limited to half the sampling frequency, is there some practically applicable way to interpolate it so that there is no danger of producing negative values? The approach should be scalable to arbitrary accuracy so basic linear interpolation won't do. By arbitrary accuracy I mean an arbitrary small difference between the band-limited continuous-time signal and the interpolation.

Perfect sinc interpolation would work because negative contributions would get cancelled by positive contributions. Ideal sinc interpolation is not possible in practice so other filters are used. This could give negative values because usually the filter impulse response has negative values. For example cubic Hermite spline piece-wise polynomial interpolation gives this impulse response:

Figure 1. Cubic Hermite spline impulse response has negative values.

Lagrange polynomials used as splines would also have negative values in the impulse response. B-splines have everywhere non-negative impulse responses, but their pass-band frequency responses deviate increasingly from the brick-wall frequency response of sinc interpolation with increasing B-spline order. The B-spline impulse response tends to a Gaussian function as the order is increased. A Gaussian function in time domain is a Gaussian function in frequency domain, so not a perfect brick-wall filter.

The obvious way to fix the problem would be to 1) take the absolute value of the interpolation output or 2) set all negative output values to zero, but I'm wondering if there is some better (inherently non-negative) method, because both 1) and 2) would create discontinuities in the derivative of the signal. This also creates spectral images.

I don't know how bad the problem stated in this question really is, because if the band-limited signal is non-negative, this also rules out some contradictory discrete signals like $[\dots, 0, 0, 0, 1, 0, 0, 0, \dots].$ The corresponding non-negative continuous signal would then have to be sinc, which is partially negative, which is the contradiction.

Shay Maymon has investigated "bandlimited square-roots" in his 2011 thesis Sampling and Quantization for Optimal Reconstruction. The idea there is to find (not easy) another band-limited signal (complex and at half the bandwidth) the squared magnitude of which equals the wanted band-limited signal. Interpolating that other signal and multiplying the result by its complex conjugate is never negative. That reminds me of this question about interpolation of the magnitude of discrete Fourier transform (DFT).


texts: Adding Salt n Pepper noise into an Image manually

texts:
When i looked in to the algorithm of Salt n Pepper noise in a Matlab i've found this,

case 'salt & pepper' % Salt & pepper noise
    b = a;
    x = rand(sizeA);
    d = find(x < p3/2);
    b(d) = 0; % Minimum value
    d = find(x >= p3/2 & x < p3);
    b(d) = 1; % Maximum (saturated) value

Which i simplified as this,

% [a, code, classIn, classChanged, p3, p4]
% p3 = density
% p4 = variance
% a = Image
% code = noise Type
% classIn = Image Class i.e unit8 , unit16, double

% -------- ALGORITHM ----------%

%     b = a;   <-- Assign b to the input image
%     x = rand(sizeA); <--- Generate random pixels from the image pixels
%     d = find(x < p3/2); <--- Find the pixels whose values are less than half of the mean value
%     b(d) = 0 <-- Implement minimum noise to them
%     d = find(x >= p3/2 & x < p3) <--- Find the pixels whose values are
%                                        greater than half of the mean
%                                        value & less than mean value
%     b(d) = 1;     <-- Implement maximum noise to them

instead of using imnoise, I am trying to add Salt n Pepper noise manually into the image using the above algorithm,

I = imread('2.jpg');
J = rgb2gray(I);
p3=0.05 %default
b = J;
x = rand(size(b));
d = find(x < p3/2);
b(d) = 0; % Minimum value
d = find(x >= p3/2 & x < p3);
b(d) = 1; % Maximum (saturated) value

I don't know where the final image is stored as when i use imshow for b it show a simple image without salt n pepper noise.
Where the final image is stored ? how can i see the it ?

P.S: What does it mean by image class ? that is unit8 , unit16 ?


texts: How would you use machine learning for peak detection?

texts:
I have a noisy signal and I'm trying to find a way to detect peaks with ML. The "peaks" are easy to find as human because they are rhythmic and have the same "general" shape but the amplitude and width of the desired peaks can vary from sample to sample.

Basically I can't just write a simple method to detect the peaks so I wanted to know if I could use ML to train a model that could learn to find the peaks more accurately. I'm new to ML so while I understand what a CNNs, RNNs and LSTM are, I don't know which would be useful in my scenario.


texts: Scaling of a rect function,drawing of signals

texts:
Im trying to understand on how the drawing of signals (time shifting,scaling...) works. This one example is confusing me, and a bit of help would be much appreciated. This is the starting Signal

We're susposed to give an $x(t)$ expression using rect, dirac and I've done so correctly;

$$x(t) = (t+2) \operatorname{rect}(t/2) + 2\delta(t-2)$$

Now we're susposed to draw a signal $f(t) = x(t/2)$;

And the rect function that is from -1 to 1 only "moves"to 2 and -2, the dirac moves to 4 that makes sense to me, but why doesn't the rect function move to 3 and -3. Thank you


texts:
I am using Python to generate a sine wave in order to cancel out part of a signal I'm analyzing. I would like the magnitude of this signal in the FFT to be equal to the magnitude of the signal that I want to get rid of. My questions are:

A. How does the amplitude of a signal translate to the amplitude of the resultant FFT?

B. Is there a way I can easily control the amplitude of an FFT using the amplitude of an input signal?

I can post my code if it helps.

This is the FFT of the signal I'm analyzing. I want the large amplitude at 1.4-ish(from the signal I generated) to be equal in magnitude to the large amplitude at 2.9-ish (the frequency I want to get rid of)


texts: frequency domain to time domain with magnitude and phase

texts:
I have a signal magnitude and phase in frequency domain. I need to have it in time domain but I really have no idea how to do it. I heard something about mirroring the signal but I'm kinda new at this. I have tried to do something like this in matlab but I'm clearly missing the point. And how do I define the time vector?
z = A .* exp(i*phase);
X = ifft(z);

I've attached the signal plot
https://i.stack.imgur.com/FfbTV.jpg


texts: Why complex signal has no imaginary spectrum

texts:
I am learning about complex sampling.

I am confused why $~e^{ j 2\pi f~ n}~$ has only a real spectrum. I would have thought the $j ~\sin(2 \pi f n)$ would produce a single spike in imaginary spectrum just like there is a single spike in real axis from $\cos(2 \pi f n)$.

I understand that the spectrum is one sided because the negative complex exponentials cancel out, but why is there not a one sided real and imaginary spectrum?

Many thanks


texts: Compress a signal by storing signal diff instead of actual samples - is there such a thing?

texts:
I am working with EMG signals sampled at 2kHz and 16 bits, and noticed that they "look smooth", that is, the signals are differentiable, and if I apply a "diff" function (numpy.diff in my case) the magnitude of the values is considerably lower than the actual samples.

So I am considering to do something like:

Split the signal into chunks of a given size;
Foreach chunk, using variable length quantity (or similar), create a byte list and:

For the first sample of the chunk, add its absolute value;
For the remaining samples of the chunk, add their difference, relative to the previous value;

This way, the smoother the signal, and the closer it is to the baseline, the more I expect to decrease the byte-size of each chunk, by decreasing the individual byte-size of a large part of the samples.

Although I suspect this would improve things for me, I also suspect that this is nothing new, and perhaps it has a proper name, and even more elegant/efficient ways to implement it.

So the question is: what is the name of this compression technique, and what are its alternatives and/or variants?


texts:
I am studying digital signal processing and our professor gave us a task of modulation and demodulation and designing an appropriate FIR low-pass filter in Matlab.

I have done the modulation and demodulation but I am stuck on the filter designing. I want to know how the coefficients of filer are determined.


texts: Number of multiplication operations needed in rational sampling rate change with IIR Butterworth filters

texts:
I am trying to get a general sense of how many multiplication operations are needed for a 4th-order Butterworth filter (IIR filter) with a rational sampling rate change. I am excluding any multiplication by 0's.
For example, for an upsampler of order 5 and downsampler of order 7 with a 99-tap FIR filter in between, the number of multiplication operations needed would be of order
$ \frac{99N}{7} $
Where $N$ is the sampling rate (per second). This is since an upsampler would effective have 4 0's for every 5 upsampled from the sampling rate of $N$.
I would like to know what the logic would be for the case of a 4th-order lowpass Butterworth filter. There are cascading designs, but there are also the canonical transpose and direct forms of implementing a Butterworth filter. What would be the difference in multiply operations (excluding multiply by 0) in each of these designs? Also, how can I calculate this? I have read somewhere essentially any IIR filter uses 6 multiplies per 2 orders of the IIR filter (4th order in this case, with $\frac{4}{2} = 2$, thus $6 * \frac{4}{2} = 12$ multiplies per second of input sample), but this is not the correct logic in my opinion.
What is the best way to think about this? I am pretty much self-teaching myself digital signal processing, so I am lost as to how to approach this.


texts: Is Savitzky-Golay filter a special case of an FIR lowpass filter?

texts: The Savitzky-Golay filter is explained here. It looks to me that any Savitzky-Golay filter can equally be considered an FIR filter. Am I right about that?


texts: Phase difference between signals sampled at different frequencies

texts:
I want to know that if it is possible to measure the relative phase difference between a signal that has been sampled at two different locations with different sampling frequencies? Also can that method be extended to undersampled cases as well?

Edit: Adding Matlab script to test possible solution (Eq.3) provided by Dan Boschen

clear all
close all
clc

Len = 768/121e6;
Fs1  = 157e6;
t1 = 0:1/(13*Fs1) :Len-1/Fs1; %Time vector for Channel 1
Fs2 = 121e6;
t2 = 0:1/(13*Fs2) :Len-1/Fs2; %Time vector for Channel 1

f=25e6; % Incoming signal frequency

phase_diff_in=0; % Modelling the actual phase difference taking In-Phase for now

% Creating signals
sign1 = cos(2*pi*f*t1);
sign2 = cos(2*pi*f*t2 + deg2rad(phase_diff_in) );
sign1 = sign1(1:13:end);
sign2 = sign2(1:13:end);

% Adding a reference cosine
sig_ref=cos(2*pi*Fs1*t2);% Fs1 sampled by Fs2
sig_ref =sig_ref(1:13:end);

% Test of phase difference formula in time domain
phi1=acos(sign1(1:256));% In first window of 256 points
phi2=acos(sign2(1:256));
phi3=acos(sig_ref(1:256));

T1=1/Fs1;
n=0:255;
phase_diff=2*pi*n*f*( ((T1*phi3(n+1))/(2*pi*n)) -T1)...
    - (phi2(n+1) - phi1(n+1));
phase_diff=wrapToPi(phase_diff);
figure;plot(rad2deg(phase_diff),'-*r')

As far as I understood the phase difference in this case should have been 0 but that is not the case. The phase difference (in deg) is as shown below:

Update: Simulating the code added by Dan

Fs1  = 157e6;
Fs2 = 121e6;
f=500e6;%25e6
samples = 400;
Len = samples;
Phi = 45;
phase_out=phase_scale(Fs1,Fs2,f,Phi,Len);
figure;
plot(phase_out)
mean(phase_out)

for the case when f=25e6 and phi=45 the following was obtained:

And for the case when f=500e6 and phi=45 the following was obtained:

The error increases significantly as the frequency is increased further.

Update #2: Simulation results after the code modifications by Dan

for the case when f=25MHz and phi=45 the following is obtained:

Which shows that the phase difference was measured very accurately.

Also for the subnyquist case as well @f=600MHz and phi=75, the following is obtained:

which shows that this works in the subnyquist cases as well. Hence the given solution works under the assumptions stated by Dan in 'Practical Limitations' section of the answer.


texts:
I implemented an interpolation method in C++ based on this equation

y(t) = x[n]*sinc((t-nT)/T)

where:
n - index of input signal,
x - input signal buffer,
t - is an actual time,
nT - time of input samples,
T - period of input signal.

I calculated it on 33 samples:1+2*16. I tested it on Matlab and C++ on VS 2010. It works great on close frequencies (eg. 24kHz to 20kHz) but when I tested situation where the input signal is 24kHz and has a peak on 8kHz and then I wanted to interpolated it to 700Hz it gave me a peak at 300Hz so it is totally wrong. I expected that the main couse of this situation is that I did not use the cutoff frequency for low pass filter in this equation above. Seriously, I do not know how to add low pass filtration to my method.

Thanks for all your suggestions or advise.


texts: Could you please confirm my method or correct me on computing the phases of my signals?

texts:
I did some measurements, and I have two received signals in two different distances.

By doing the fft on the signals, I have the spectrum, then by using the angle() function in MATLAB on my fft, I computed the phase for each frequency.

I understand that the angle() function just give me the angle of each point of fft. I'm confused, is this the phase of my signal? Is this the phase relative to what signal?

As far as I know, phase of a signal, a Line-of-sight signal, after traveling a distance of L is computed as this:

-2*pi* f*(L/c)

c: is the speed of light.

The number that the angle(fft()) gives me is very different. Not even close to the above formula.

And I did all of this in MATLAB.


texts:
As a project requirement, I am supposed to modulate a complex baseband signal to a carrier frequency of 2.4 GHz. These are the steps I followed for the baseband modulation and demodulation:

Generate Data -> OQPSK Modulation -> Up Sample -> Pulse shaping -> AWGN Channel -> Matched Filter -> Down Sample -> OQPSK Demodulation

The sampling frequency at baseband is 12 MHz. Now if I were to modulate this baseband signal to a carrier at 2.4 GHz, I believe that these are the steps to be followed:

%Generate Carrier
t = (0 : 1/Fs : (length(baseband_Signal)-1)*1/Fs)
carrier = sqrt(2)*exp(1i*2*pi*Fc*t);

where sampling frequency $Fs=12MHz$ and carrier frequency $Fc=2.4GHz$

% Frequency upconvert to passband.
xUp = real(baseband_Signal.*carrier);

% Downconvert to baseband
yACI = rx_sgl.*conj(carrier);

When I perform upconversion, I expect to have 2 peaks centered at +2.4GHz and -2.4 GHz when I plot the frequency spectrum. However I see that $xUp$ is centered at 0.

I am not sure about what I am missing. Should the complex baseband signal again be upsampled ?


texts: Reading about singular value decomposition (SVD) in the context of signal processing applications, one can separate the signal from the noise into orthogonal subspaces. On the surface this sounds like using the decomposed output of the SVD is more optimal than matched filter (or correlation) for detection, demodulation, timing estimation, etc. Is this true or is there more nuance to this problem?


texts: simulink discrete integrator block problem

texts:
I am testing the Simulink discrete integrator block as shown in the picture below. The simulation is fixed step with $T_s = 1e-5$ and the cosine wave block setting and the scope display is also shown below.

For the discrete integrator block, the gain value K is set to be 1 and the integration method was chosen to be "Accumulation: Forward Euler". However, how come the magnitude of the result is so big as shown in the top part of the scope picture? I thought integrating a cosine wave will produce a sine wave with the same magnitude. If the integration method is chosen to be "integrator: Forward Euler", the result becomes very small, as shown below. It just can't produce a sine wave with magnitude of 1. How to get the correct result?


texts: Recoverable continuous function after average sampling

texts:
I have a personal interest in this question, but I don't know if it is unfounded.
Suppose I have a continuous function x (t), which is sampled using an average sliding window, which performs a continuous uniform integration. The start of the scan is not fixed, that is, the signal could be scanned at any time. Is it possible to generate a function that is resistant to this procedure and can be recovered with an acceptable error from the samples obtained? What are the limitations in the generation of this function?
Thank you!


texts:
Given a signal $$\sin\left(\frac{2\pi}{\sin(2\pi/Tx)x}\right)+n$$ (where n is noise) is there a way to isolate and decompose the number, frequency, and causal (carrier-modulator) order of the oscillators (here 2)?

(Similar to How do I find amplitude modulation effects in natural signals but more formal).


texts: Use of spectrograms in using deep learning on EEG related problems

texts: I am currently working on a classification problem related to EEG signals. My professor asked me to convert those signals into their spectrograms and try various convolutional neural networks on them. However, he didn't provide any reason or justification for his suggestion. What could be the benefit of using spectrograms for such problems instead of going for something like RNN or LSTM since we are dealing with sequential data here?


texts:
I was looking for a more efficient way of finding the magnitude and phase of a signal at a certain frequency without performing an FFT because it produces more information than I need and I came across this comment on the Matlab message board:

http://www.mathworks.com/matlabcentral/newsreader/view_original/250421.

If I am sampling at 1500Hz and my signal is 100 points and I am looking for the phase and magnitude at 15Hz, would this be suitable method for obtaining the information I need? Also, why does the OP suggest windowing the signal as well?


texts: How to remove wind noise from video file

texts:
I've some video footage where wind noise is drowning out background noise I want to be clear.

Is it possible to do some DSP to remove/attenuate this wind noise in order to hear the engine noise?

I am not well versed in DSP, but my initial assumption would be;

Extract audio from video file
Calculate FFT of audio file so as to analyse the frequencies at which wind is present
Design and implement a filter which removes unwanted noise from audio file
Modify original movie file by replacing the existing audio with the "cleaned" audio

Any help appreciated.


texts: Why does over-modelling an adaptive AR NLMS filter fix sharp spikes?

texts:
I just simulated an auto-regressive second-order model fueled by white noise and estimated the parameters with normalized least-mean-square filters of orders 1-4.

As the first-order filter under-models the system, of course the estimations are weird. The second-order filter finds good estimates, although it has a couple of sharp jumps. This is to be expected from the nature of NLMS filters.

What confuses me is the third- and fourth-order filters. They seem to eliminate the sharp jumps, as seen in the figure below. I can't see what they would add, as the second-order filter is enough to model the system. The redundant parameters hover around $0$ anyway.

Could someone explain this phenomenon for me, qualitatively? What causes it, and is it desirable?

I used step size $\mu=0.01$, $10^4$ samples, and the AR model $x(t)=e(t)-0.9x(t-1)-0.2x(t-2)$ where $e(t)$ is white noise with variance 1.

The MATLAB code, for reference:

% ar_nlms.m
function th=ar_nlms(y,order,mu)
N=length(y);
th=zeros(order,N); % estimated parameters
for t=na+1:N
    phi = -y( t-1:-1:t-na, : );
    residue = phi*( y(t)-phi'*th(:,t-1) );
    th(:,t) = th(:,t-1) + (mu/(phi'*phi+eps)) * residue;
end

% main.m
y = filter( [1], [1 0.9 0.2], randn(1,10000) )';
plot( ar_nlms( y, 2, 0.01 )' );


texts: How should buffer value at N/2 get initialized prior to performing an inverse FFT on a buffer, size N, that contains a complex frequency response in only values 0 to N/2-1?  It is pretty clear the buffer needs to be made symmetrical for prior to calculating the inverse FFT per table 12-1 of http://www.dspguide.com/ch12/1.htm.  What is not clear is how to initialize the value at N/2.  Is it related to value at 0?


texts:
When rectifying a signal by squaring it, is it necessary to upsample by a factor of 2?

Why is this?


texts:
I've stacked on solving this problem and can't find an answer in the web,  so I decided to ask a question here.
There are the problem:
Preselector of boadcast receiver contains of one oscillation circuit, that has a bandwidth 8 kHz at LW range. Quality factor at SW range is 120. How will the selectivity of the receiver on the image frequency channel be changed if we will switch from LW range($f_{RF}=280 kHz$) to SW range($f_{RF}=12.04 MHz$)?
The answer is: It will decrease for $17,4 dB$.
The selectivity os given by formula $\sigma=10\lg(1+\xi^2)$, where $\xi=Q(\frac{f}{f_0}-\frac{f_0}{f})\approx\frac{2Q\Delta f}{f_0}$.
And for the image frequency channel $\Delta f=2f_{IF}$
I'm understand what is going on and how to solve this problem if the $f_{IF}$ if done, because I can get the quality factor on LW range from the bandwidth, but there are no $f_{IF}$ in the problem.
Any ideas will be appreciated!


texts: Peak Signal to Noise Ratio (PSNR) in Python for an Image

texts: How do I calculate the maximum signal to noise ratio (PSNR) in Python? Is there any library that can calculate PSNR for Image?


texts: Faster way to find the quantization level of a sample

texts:
Given a sequence of real numbers that represent a sound signal i.e [1.8, 2.2, 2.2, 1.9, -1.5, -0.7], we must quantize this sequence by dividing the range [-4, 4] to 32 equal parts.

If we had to do this with the traditional way we would need to write 32 different ranges

[-4   , -3.75]

[-3.75, 3.5]

.
.
.

and check in which one the real value is in.

Is there an another way to achieve this? Maybe with a normalization function or what?


texts: how to guess what the interference source it could be

texts:
At work, we are in an industrial area, and there is always a signal that I receive inside my office which interferences with my tests.

This signal is in the 433 MHz and when I see the spectrum, there are 2 peaks around this center frequency which are shown sort of constantly pulsing, like if someone were pressing a button every 5 seconds (but this is a 24h interference, so it must be something else).

Also, because of the strength of this signal, I wonder if it's something inside my office but we only have routers, wireless phones...I can't think of any wireless device inside my office that could be transmitting at this frequency.

I only can thing about the garage door openers but, as this interference is constantly there, it must be something else.

What sort of device could be doing this?


texts: complex multiplier in divide and combine FFT

texts:
I am studying radix 2 algorith from Proakis' book.

But I'm a bit confusied why 1st DFT $G_1$ is not multiplied by complex entity while 2nd DFT $G_2$ is being multiplied by complex entity $W$ as shown highlighted in attached figure.


texts: Why is particular solution zero for an impulse excitation signal?

texts:
We were being taught the impulse response for a series RC Circuit- consisting of simply one resistance, one capacitor, and an impulse excitation all in series.
I get that the homogeneous part of the solution will be of the form:
$$x(t) = Ke^{\frac{-t}{RC}}$$
And we can find the value of $K$ using the initial conditions.
However, when it came to the particular solution part, we were straightaway told that since the excitation is an impulse, the particular solution will be equal to $0$. I'm not able to get why is this the case.
Also, I am not getting exactly how to calculate the particular solution for a given network, say RC network.


texts: Python MNE - reading EEG data from array

texts:
I have EEG data that comes in the form of a 3D numpy array (epoch * channel * timepoint). timepoint is a 256 element array containing each sampled timepoint (1s total, at 256Hz). epoch is an experimental trial.

I'm trying to import the numpy array into a form Python-MNE (http://martinos.org/mne/stable/mne-python.html) understands, but I'm having some trouble

First, I'm not sure if I should be importing this raw data as a RawArray or an EpochsArray. I tried the latter with this:

ch_names = list containing my 64 eeg channel names
allData = 3d numpy array as described above

info = mne.create_info(ch_names, 256, ch_types='eeg')

event_id = 1

#I got this from a tutorial but really unsure what it does and I think this may be the problem
events = np.array([200, event_id])  #I got this from a tutorial but really unsure what it does and I think this may be the problem

raw = mne.EpochsArray(allData, info, events=events)

picks = mne.pick_types(info, meg=False, eeg=True, misc=False)

raw.plot(picks=picks, show=True, block=True)

When I run this I get an index error: "too many indices for array"

Ultimately I want to do some STFT and CSP analysis on the data, but right now I'm in need of some help with the initial restructuring and importing into MNE.

Whats the correct way to import this numpy data that would make it easiest to complete my intended analyses?


texts: I have a two data set files data1 and data2 each containing 5000 points having 3 columns representing x,y(coordinate (centroid)),angular movement(theta) resepectively for two robot movements obtained by numerical programming and solution of a circuit representing robot equations.I need to simulate a graph for coverage with obstacle and without obstacle in a grid like fashion and then calculate the percentage of coverage where A is the coverage for each cell, A(i)=1 if cell i is covered otherwise 0. I went through this link SO but the answer is too obscure and I cannot follow. Can somebody please tell me how to generate these two graphs?Thank you


texts: Confusion with regards to STFT and MFCCs

texts:
I'm wondering if someone can clear the following up for me as I'm having a difficult time understanding the concept.

I'm looking to implement and learn about the processes and uses using short-time Fourier transforms. I have had experience using DFT's and wanted to expand my knowledge.

From research that I have read online, as well as in some text books I understand that the STFT is calculated by the following:

Take the inputted signal of $x[n]$
Multiply $x[n]$ by a window function (e.g. Hamming, Hanning, Triangular, .., Barlett) by sliding the block size over the time axis
Compute the DFT for the resulting values of step 2

Is this correct? I know the following steps above are for usually a continuous-time STFT and in the case of a Discrete-time STFT the data can be broken up into chunks.

The other thing that is confusing me:

For a project last year, I had to implement a MFCC algorithm. Inside this algorithm I computed the "Mel Scale" in Triangular filter bank and then multiplied each of these values against the FFT (DFT) computed on the inputted signal in Discrete-Time. Is this therefore a use of the STFT, just using the Mel-Frame over Hamming etc.. Or, would I still need to multiply the result of the FFT (DFT) with the Triangular filter bank if I was to use a STFT in order to calculate the MFCC values?

Thank yo for reading.


texts: Best way to implement variable phase shift on FPGA?

texts:
I have an incoming digitally sampled sine wave pulse, so my FPGA has the ADC level from an I and Q channel. I want to be able to shift the phase by some arbitrary amount.

What's the most efficient way to do this?

The most obvious way to shift phase would be just to multiply by $e^{i\phi}$, but that would be slow and wasteful. I've also seen someone recommend using an allpass filter with some phase shift.

Is there any other options I should be considering? Or is the allpass filter the best way to do this?


texts: What are the types of noise (FM,AM, MSK modulated) that I can use for jamming GSM signals?

texts: I need to jam GSM signals in-order to prevent mobiles in near-by vicinity from connecting to the BTS. I am employing a noise based jamming strategy. What type of noise will be most appropriate for this task?


texts:
I'm conducting experiments on several n-ary Psk signals of differing symbol rates. I'm comparing detection levels of different signal detection approaches on sets of PSK signals. I'm using a cosine-roll-off filter with a roll-off-factor of between 0.2 to 0.5.

To determine the appropriate signal to noise ratios, I had to decide which bandwidth I was gonna use. For reasons of practicality, I'd like to use the approximation symbolrate[bd] equals bandwidth[Hz]. It was proposed in chapter 6.2.4 of the book:

Communications and Information Systems; Michael John Ryan,Michael Frater

Do I act reasonably and doesn't it skew my evaluation with inappropriate error levels?


texts: I already have the value of Image Rejection Ratio (IRR) in dB for IQ imbalanced modulator and I want to simulate this effect using MatLab. How can I obtain $\alpha$ and $\theta$ to calculate the imbalanced gain for I & Q branches ??


texts: Spectral estimation in Spectrum analyzer

texts:
I understand that there are various techniques for spectral estimation of the signals.
Methods such as parameteric(AR,ARMA etc.), non-parameteric( periodogram , welch etc..).
But I fail to understand which one is used by spectrum analyzer.
I am unable to relate the theoretical methods with spectrum analyzer.
Can anyone explain please?


texts: How the Saturated Pixels violate the assumption that the image formation process in linear?

texts:
I am working in digital image restoration field and i have studied this is many research papers that presence of saturated pixels violates the assumption which is initially assumed by many blind deconvolution algorithms, that the image formation process is linear...

Can anyone please explain that how this rule is violated.?


texts: Remove a frequency and all its multiple (harmonics)

texts:
In a soundfile, I'd like to remove a frequency and all its harmonics. What is the best solution in your opinion ?

Do a loop with 20 notch filters, one after another.
// f is the fundamental frequency
for (i=0; i++; i<21)
{
file.notchfilter(f*i);
}

Question : does each notch filter create some unwanted side-effects, and the fact that we use 20 filters one after another will create very high unwanted artefacts?
Is it possible to design a filter that will remove 20 frequencies in ONE STEP ?
Then we wouldn't need to do a loop, just one pass wil be okay?
Is it called a Comb filter ?
Other method ?


texts: Cross power spectral matrix from PSDs and coherence using matrix notation

texts:
I have a vector containing all the PSD of a series of signals

$\underline{s} = \left[ |S_{x_1x_1}| |S_{x_2x_2}| |S_{x_3x_3}| ...\right]^T$

I also have a matrix $\gamma^2$ containing all the coherences

$\gamma_{ij} = \frac{|S_{x_ix_j}|^2}{S_{x_ix_i} S_{x_jx_j}}$

I would like to obtain the matrix of the cross-spectral densities, i.e. a matrix $\mathbf{S}$ so that
 $S_{ij} = |S_{x_ix_j}|^2$

Is it possibile to obtain it from $\underline{s}$ and $\mathbf{\gamma}$ using matrix notation?

The only way I found at the moment is:

$\mathbf{S} = (s s^T) \circ \gamma$

where $\circ$ indicates the Hadamard (element-wise) power and multiplication. This is however seldom used and usually element-wise multiplication is avoided in linear algebra.

Is there another way to do it?


texts: Cross spectrum using `cpsd` (matlab) of 2D matrices

texts:
I'm struggling with cpsd for 2D matrices. In my test, I've generated two 100x100 matrices (sine waves) s1 and p1 with a small phase difference between them. I'm trying to identify the propagation direction using this phase difference, hence using cpsd

[Psp F] = cpsd(s1,p1....).

My guess is that cpsd is the one for the job. However I'm not sure if I'm doing the right thing. When the crests of s1 and p1 are orthogonal to the axis, my script works just fine. However, if they have some angle, instead of two components, let's say in the 1st and 3rd quadrant, I end up with 4 components, in all quadrants. Not good...

I suppose there's something wrong with the way I'm using cpsd for matrices. Anyone out there could give me a hand? Many thanks.

This is my code so far

[numrow numcol] = size(s1);
for c = 1:numcol
    [Psp1(:,c),F]= cpsd(s1(:,c),p1(:,c),[],[],[],dt,'twosided');
end
for l = 1:numrow
    [Psp2(l,:),F]= cpsd(s1(l,:),p1(l,:),[],[],[],dt,'twosided');
end
Psp=Psp1*Psp2;
F1=F-1/(2*dt); %
[XF YF] = meshgrid(F, F);
[XF1 YF1] = meshgrid(F1, F1);


texts: Gain at critical frequency for Butterworth and Chebishev filters

texts: In Chebyshev filter the critical freq occurs with gain Gp But in butterworth filter it is always at -3dB; why does this happen? Can anyone explain intuitively on this?


texts: What are some reliable sources which I can use to validate my code for simulating MPC?

texts:
What I have done
I have written a code in Mathematica from scratch for simulating quadratic cost standard MPC problems for constant reference tracking.
My attempts
YALMIP's website has a solved standard MPC example which should be reliable but it has a linear cost, which requires quite different functions and a different code (some parts are same as in quadratic case but some are different) to simulate and thus it does not validate my original code for quadratic costs.
I have tried to validate my code against a number of solved examples in research papers and lecture notes but the outcome is a mixed bag. For some my results match perfectly, for some the results for piece-wise constant reference tracking match only for some constant references (and an attempt to track the other references yields the error that there are no points that satisfy the constraints) and for some I get the error that there are no points that satisfy the constraints.
The problem
I am in a fix here. For each case above, I have tried at least two sources. I have found an outright mistake in one of the sources, a research paper, and this makes it difficult to know if any other source can be relied on or not.


texts:
I have some images of a background, and of steam that I am subtracting to extract just the steam. Unfortunately, because the camera conditions are not perfect the differenced image is always a few pixels off generating sharp lines where features in the background don't match up pixel-by-pixel like so

I'm interested in segmenting out the steam (foreground noise), from the background. I think if I can automatically remove the continuous lines from this image it would be sufficient for my purposes.

Unfortunately it's been almost ten years since my image analysis work and I don't remember what you do to detect continuous lines. What are some simple well preforming algorithms to do this?


texts: Fourier Transform Signals - Time Transformations

texts:
I was going over some review problems and came across an interesting one.

Using the techniques of (linearity, time shifting, and time scaling) what are some approaches I could use to turn the fourier transform from the left graph to the right?

X(ω) = 4* (sin w)^2 / w^2

My understanding so far:

Using the time shift property for fourier transforms I know the we can use
x(t -1) to move the graph to start at t=1
We can use x(t/3) to extend or "stretch" the base to the to t=7

What is eluding me is how to form that "middle" or clipped top of the z(t) graph.
Any insight?


texts: What is an intuitive explanation of cepstrum for a 12-year-old?

texts:
In Feynman's spirit of "If you truly understand a concept, you should be able to explain it clearly" and the explanation should make solid sense.

There are two main challenges:

The axis of "Quefrency" is not a straightforward/natural concept, how to describe and make sense of what it is?
Why is the cepstral view an efficieny way of "elevating" the main information contained in the signal? We can see that the vector is more succinct, but where does this efficiency come from?


texts: RC high pass filter differential equation

texts:
I have the question to derive:
$$RC \frac{dy(t)}{dt} + y(t) = RC \frac{dx(t)}{dt}$$
which describes a high pass filter or $CR$ circuit with $x(t)$ the input and $y(t)$ the output.
I am really struggling to derive this
I start with Kirchhoff's voltage law
$$x(t) = Ri(t)+y(t)$$
Use the differential equation for capacitor current. (in series so resistor and capacitor current are the same)
$$i(t)=C \frac{dy(t)}{dt}$$
I then substitute this in to Kirchhoff's voltage law
$$x(t)=RC \frac{dy(t)}{dt} +y(t)$$
This is the same workings out as for a low pass filter (RC circuit). Is this correct?
If so how do I get to the initial equation?


texts:
I'm dealing with a brain stimulation.

My EEG contains large spikes (>3 mV) when stimulation occurs. I want to remove the 50 Hz noise, but I can't use notch filter because it creates large artifacts.

Do you think that Kalman Filter could remove the 50 Hz noise?

Thanks for your comments. This is my registration

If i apply a notch filter i will have ringing artifacts due to the large spikes.If so I can't use the signal anymore. So my point is : i have to remove the 50hz component without using a notch filter.( I cannot delete the spikes with interpolation) My idea is to estimate the signal rappresenting the 50 hz component  with the kalman filter and then subtracting from the original signal


texts:
Basic questions: What's the "correct" way to deconvolve two causal FIRs in the frequency domain (i.e. using the FFT), neither of which may be minimum phase but both may be considered to have stable exact inverses? By "correct", I mean, how does one avoid time-aliasing errors? If there is no "correct" way, what is the best way?

Some context: Given two causal FIRs, it is easy to avoid time-aliasing when convolving them in the frequency domain by first appropriately padding the IRs with trailing zeros (i.e. performing linear convolution via the FFT as discussed in standard DSP texts). However, it's not clear to me how to do something similar for deconvolution when the FIRs may not be minimum phase.

Deconvolving a minimum phase FIR out of a causal FIR (minimum phase or not) seems straightforward since we can compute the exact inverse of the minimum phase FIR (which we know will be causal and stable), and then linearly convolve this inverse with the other FIR.

However, since the exact inverse of a causal FIR that's not minimum phase is non-causal (assuming the inverse is BIBO stable - we force the inverse to have a finite length), performing linear convolution using the FFT after computing the inverse does not make sense since padding a non-causal inverse with trailing zeros doesn't make sense.

Related example with somewhat baffling results: In the figures below, I linearly convolve $x1$ with $x2i~(i = 1,2,3)$ to produce $yi$. Then, I deconvolve $x1$ out of each $yi$ to try and reproduce each $x2i$. To do this, I pad $x1$ with trailing zeros until the length matches that of $yi$, compute the exact inverse of the padded $x1$, and then circularly convolve the inverse with $yi$. The relevant difference between the plots across columns are the signals $x2i$ which have a different number of trailing zeros ($i$ denotes the number of trailing zeros). Can anyone provide some explanation as to why the results seem reasonable in columns 1 and 3 but not in column 2? Also, this might seem to be a reasonable way to deconvolve two signals and yet the results obtained don't always appear satisfactory (as this example shows). Any thoughts?

Note: In the plot titles, I have used $*$ to denote linear convolution and $/$ to denote deconvolution (all done in the frequency domain as described above).

Thanks!


texts: What is this effect/filter in given images?

texts:
Here are two images. What kind of filter/processing is done here?
It looks like an emboss but the bumps are smooth.

1) Original :

2) Processed :


texts: Are there any modulation schemes (or other situations) where sampling at or very near twice the carrier frequency would be better than IQ or quadrature sampling using a local oscillator at or near the carrier frequency for SDR demodulation?


texts: Impulse response of forward difference cascaded with one sample delay

texts:
Below is the excerpt from Discrete Time Signal Processing by Alan Oppenheim.

I don't get how $(\delta[n+1] - \delta[n]) * \delta[n-1])$ becomes $\delta[n] - \delta[n-1]$. The convolution sum operator "*" here is defined as below.
$$y[n] = \sum_{k=-\infty}^{\infty} x[k]h[n-k] $$
How do I substitute the definitions of forward difference and one sample delay in this and arrive at the impulse response the author has given?


texts: What does it mean that one frequency aliases another?

texts:
So I was watching a video on sampling and the professor asked which frequencies will alias to 0 and which to $B$ if we are sampling at frequency $f_{s1}$. What does mean that two frequencies alias each other ?

The answer is that $0, f_{s1},2f_{s1},3f_{s1},\ldots, nf_{s1}$ will alias to 0 and $f_{s1}+B, 2f_{s1}+B, \ldots, nf_{s1} + B$ will alias to $B$.

Could anyone explain this please? The image shows the frequency spectrum in question.


texts:
In my signal processing class I have some basic tasks working with sinus signals in Matlab. I have two sinusoidal signals with 1kHz and 3kHz and sampling rate of 16kHz. The next task says 'Create a linear combination of the two signals and save the result in a new vector'.

I don't really understand what's meant by the linear combination. I remember something vague about having to add them with weights but I'm don't quite comprehend what that does or how to do it.


texts: How to perform downconversion using a DSP

texts:
I read that to perform downconversion, first, I need to generate a sine wave with the same phase and frequency as the transmitter's and then multiply the signal received by the generated sine wave. I don't know how to do that, I think I have a way to generate the synchronized sine wave, but  I don't know how to perform the multiplication.

I read that I can use a mixer or a nonlinear circuit to perform the multiplication, but I have a DSP and it doesn't have a frequency mixer available. Can I do it digitally ? The carrier has 125kHz and the modulating signal has 20kHz. I can use a sampling rate up to 1MHz and the DSP can reach up to 70 MIPS.

I'm sorry if my explanation is incomplete, I am really studyng the problem, so if you need more information you can ask and I will try to make the question more complete.

With my best regards,

Daniel.


texts: Is Cross correlation between stationary random signal and noise zero for all lags?

texts:
If $x(n)$ is a stationary random process (let say small 20ms window of speech which can be assumed stationary) and $v(n)$ is said to be uncorrelated noise, can we say that cross correlation of $x(n)$ and $v(n)$ is always zero, for all lags? E.g.

$$r_{xv}(\tau) = E\{x(n) v(n+\tau)\} = 0, \forall \tau$$

Or,

$$\sum_{n=0}^{N}x(n) v(n+\tau) = 0, \forall \tau$$

where N is the signal length?


texts: Straight lines with an edge detection algorithm

texts: I've been thinking of an algorithm for object recognition but it would relies heavily on straight non-noisy lines and as far as I know this is hard to obtain with edge detection algorithms. What's the best edge detection algorithm to obtain lines as straight as possible?


texts: Can a wavelet transform give time dependent phase of sinusoids in signal

texts: I have a signal which contains sinusoidal components that oscillate at different frequencies. I think the phase of the sinusoids is changing with time. I could do a Fourier transform on small chunks of the signal and the phase of the sinusoids from that, but is there a better way to do this with a wavelet transform? I'm new to the concept of wavelets. From what I understand, it's good at doing what an FFT does but with time resolution. Could you get the phase of a sinusoid as a function of time from a wavelet transform ?


texts: Kalman filter for position and velocity: introducing speed estimates

texts:
Thanks to everyone who posted comments/answers to my query yesterday (Implementing a Kalman filter for position, velocity, acceleration ).  I've been looking at what was recommended, and in particular at both (a) the wikipedia example on one dimensional position and velocity and also another website that considers a similar thing.

Update 26-Apr-2013: the original question here contained some errors, related to the fact that I hadn't properly understood the the wikipedia example on one dimensional position and velocity.  With my improved understanding of what's going on, I've now redrafted the question and focused it more tightly.

Both examples that I refer to in the introductory paragraph above assume that it's only position that's measured.  However, neither example has any kind of calculation $(x_k-x_{k-1})/dt$ for speed.  For example, the Wikipedia example specifies the ${\bf H}$ matrix as ${\bf H} = [1\ \ \  0]$, which means that only position is input.  Focussing on the Wikipedia example, the state vector ${\bf x}_k$ of the Kalman filter contains position $x_k$and speed $\dot{x}_{k}$, i.e.

$$
\begin{align*}
\mathbf{x}_{k}  & =\left(\begin{array}[c]{c}x_{k}\\
\dot{x}_{k}\end{array}
\right)
\end{align*}
$$

Suppose the measurement of position at time $k$ is $\hat{x}_k$.  Then if the position and speed at time $k-1$ were $x_{k-1}$ and $\dot{x}_{k-1}$, and if $a$ is a constant acceleration that applies in the time interval $k-1$ to $k$, from the measurement of $\hat{x}$ it's possible to deduce a value for $a$ using the formula

$$
\hat{x}_k = x_{k-1} + \dot{x}_{k-1} dt + \frac{1}{2} a dt^2
$$

This implies that at time $k$, a measurement $\hat{\dot{x}}_k$ of the speed is given by

$$
\hat{\dot{x}}_k = \dot{x}_{k-1} + a dt = 2 \frac{\hat{x}_k - {x}_{k-1}}{dt} - \dot{x}_{k-1}
$$

All the quantities on the right hand side of that equation (i.e. $\hat{x}_k$, $x_{k-1}$ and $\dot{x}_{k-1}$) are normally distributed random variables with known means and standard deviations, so the $\bf R$ matrix for the measurement vector

$$
\begin{align*}
\mathbf{\hat{x}}_{k}  & =\left(\begin{array}[c]{c}\hat{x}_{k}\\
\hat{\dot{x}}_{k}\end{array}
\right)
\end{align*}
$$

can be calculated.  Is this a valid way of introducing speed estimates into the process?


texts:
Say we have a function of time ($S(t)$) of the length $T$, and then a customized impulse response (say $I(t)$) of the length $T+N$. The question is, when $S(t)$ is convolved with $I(t$), what are the possible frequency domain (or alternately time domain) representation of all of the possible resulting functions?

For the case when the convolution is performed circularly, the solution appears immediately clear. If $S(t)$ has zero magnitude components in the frequency domain, they will stay zero no matter how the impulse response is configured. Thus, it's impossible to for example, convert a sinusoidal waveform into any non-sinusoidal waveform using circular convolution.

But the result is a lot less clear for the case when linear convolution is performed. If we limit our horizon to the period of $T$, and completely disregard the tail that the impulse response will generate (for cases where $T>1$, and length of the impulse response $ > 1$ ), what will the answer be then? Are there any limits, or can you transform any function into any other function (provided the tail is thrown away), using linear convolution?

(EDIT: for non-casual impulse responses the "pre-echo" is also thrown away, time shifting is allowed also for such cases).


texts: Relation between time domain, DTFT domain and frequency domain

texts:
Problem

The sampling frequency of a continuous-time signal is $S$ kHz, what does $\frac{\pi}{4}$ radians/sample in DTFT domain represent in Hz in frequency domain? Prove the relationship.

Doubts

I am confused with DTFT domain and frequency domain and I do not how to put all these domains together.


texts:
I am struggling with understanding the phase contribution of each individual pole. Let's say we have a system (minimum-phase system if it makes a difference) and it has poles located at:

and

What is the analytical relation between the phase contribution of these two complex poles?


texts: What is state of the art in multivariate 2D image display?

texts:
If I have a 2D scalar image, it's trivial to display it with MATLAB or matplotlib type tools with the pixel values mapped via a color map, and there is any amount of literature (one example, another one) describing the perils of those tools' default "rainbow" colour-maps and the advantages of various other schemes (color ramps constructed to be more suitable for the significant proportion of the population which is colorblind, for example).

However, if I have two 2D scalar images and want to fuse/overlay them in some useful way (rather than just displaying them side by side), it's less clear what the best approach is.  The naive approach of using one channel as red and one channel as green in a combined image kind of works, but is pretty ugly looking.  Mapping one channel to luminance and one to hue in an HSV color space also kind of works, but it becomes impossible to identify the hue of the dimmest pixels (unless luminance is kept quite high, but that destroys contrast).

So my question is: is there some recognized and generally considered good/best way of visualizing such fused image pairs?  Or does it all depend on the details of the data and the specifics of the "visualization problem" to be solved?

NB I'm not considering the case of visualizing images of vector quantities or complex numbers; this is more about independent but spatially registered scalar image data from different sensors/"modalities" or GIS/mapping systems where the viewer wants to get some insight into spatial correlations between the imaged quantities.


texts:
I am new to DSP, and I'm self studying. I am confused about the magnitude of $e^{j\omega}$ - where $\omega$ is the normalized angular frequency - when we are on the unit circle.

According to the text book, it's $|e^{j\omega}|=1$. But why is that? I'd really appreciate it is someone could please kindly explain this to me.


texts: How to downscale an image with interpolation by a factor larger than interpolation kernel size?

texts:
I'm experimenting with image resizing techniques and algorithms. Specifically, I'm significantly downsizing images, e. g. from 2048x1536 to 64x48 - 32 times.
Now, say I'm using a 4x4 kernel. Right now I'm upscaling the image to be by a factor of power of 4 larger than the target (2048x1536 becomes 4096x3072), and then process the image 3 times, each time obtaining a new image 4 times smaller than the source. Seems ineffecient. Is there a trick to it?


texts:
I measure vibration along a beam, 74m long. I have 2 sensors either side of the beam and measure vibration for 1 second ($t= 1s$) at 9000 Hz sampling frequency. I measured speed to 12000 m.s.

I'm picking up features along the beam, representing discontinuities along the beam. I would like to convert the time vector to distance, so that everything which happens in $t$ is now distance along the beam ($x$ axis) and $y$ axis represents the discontinuity feature.

This is in MATLAB.


texts: Finding the transfer function of a discrete signal described by two equations

texts:
A discrete time system is described by the following system of equations.

$$q[n] = \big(x[n]-\frac k4q[n-1]\big)$$
$$y[n] = \big(q[n]-\frac k3q[n-1]\big)$$

Find the systen function and then find the values of $k$ for which the system is stable. Also find the values of $k$ for which the system is of mimimum-phase.

Can anybody give me some guidelines as to how to handle this system of equations in order to find the system function?

This is an unsolved exercise given by my professor.


texts: Comparing the FFT to numerical integration in Matlab

texts:
The result of calculating the fourier transform using numerical integration is:

the result of using Matlab's FFT is:

So where did I go wrong here? I know the FT of a Gaussian should be another Gaussian.  Shouldn't they be the same? Furthermore I would have thought the
Edit:
Following a suggestion using the fftshift function gets me to the following plot, but this looks very different to the original fft despite documentation stating that it just centres the result around zero frequency.

which is more resembling what I would expect.
omega = linspace(-30,30,1000);

t = linspace(-10,10,1000);
y = exp(-t.^2./2)/sqrt(2*pi);
Yft = fft(y);

% Pre-allocate results of integration.
Y = zeros(size(omega));

%Integrate.
for k=1:length(omega)
    Y(k) = trapz(t, exp(-t.^2./2)/sqrt(2*pi).*exp(1j*omega(k).*t));
end

figure()
plot(omega, real(Y), omega, imag(Y));

figure()
plot(omega, real(Yft), omega, imag(Yft));


texts:
I'am pretty new to this topic so i dont have much experience with DSP.
I want to filter(highpass) a WAVE file, i have already programmed the FFT, and now want to filter the FFT vector via convolution (which is a multiplication in frequency domain ?!).

1) Can anyone give me an easy equation for such a filter?

2) If i have an equation i compute it, also Fourier-Transform it, an multiply it with my FFT vector and inverse Transform it an i have applied the filter?!


texts: Proof of time-invariance of continuous-time system

texts:
I have a system with the following input/output relation:
$$ y(t)=x(-t) $$
and I want to prove (not graphically/draw) that its not time invariant (TI).
I tried to write down $y(t-T)$ and compare it to the response to $x(t-T)$ but I get both terms equal which shouldn't be the case.


texts: How to objectively assess the quality of a medical image(CT,PET,MRI)? What factors should be considered?


texts:
I have done my receiver design with QPSK modulation. I am trying to communication between two systems at L band. My data rate is 1 Mbps. As my data rate is high I have allocated 64 bit for timing recovery.

Timing recovery has farrow interpolation filter followed by Gardner timing recovery, loop filter (2nd order) and interpolation control for the farrow interpolation filter. I am doing 2 samples/symbol analysis.

What I am assuming is the fractional value to the interpolation filter will stabilize once 64 bit allocated for timing recovery is over. But what I observe is the fractional value varies even after the 64 bit allocated is over because of which my data sometimes is wrong.

can the fractional input to the interpolation filter vary?
if yes why is the data sometime is wrong even if it is just direct board
to board communication?


texts:
If I have a noisy signal which I know is approaching the saturation limit of my sensor, then the measurement will sometimes be saturated. It seems that if I try to estimate the true data by averaging the noisy measurement, then the mean will be systematically shifted from the true value (assuming zero mean white noise).

In the example below, when the signal approaches the saturation limit of 100, the moving average is systematically below the true signal. The moving median would solve this particular problem nicely, but I would like to avoid using a median filter for other reasons*.

Is there any "standard" or well known way of dealing with this kind of systematic error which can be expected in advance? (Apart from median filtering.)

I'm interested in the general case, but in case it helps to clarify things, I will include my specific application. I am trying to categorise the brightness of various regions in an image into discrete levels. For black regions, saturation is at 0. Due to noise, most pixels are >0, none are <0, so the average is >0, which sometimes fools my routine. Over the rest of the dynamic range, I can differentiate quite well between N and N+1 levels of brightness. Differentiating between the 1st and 2nd level seems to be more challenging, because of the saturation issue. In the dark (not black) regions, the noise is symmetrically distributed, so I guess I could test for skewness? But I think this would not be very robust for a small number of pixels.

*the other reasons being that I'm actually using a weighted average or Gaussian filter to take advantage of spatial structure within the image. I am not aware of a "weighted median" filter, and I'm assuming at this stage that such a thing doesn't make a whole lot of sense.


texts: Why is the difference of gaussians scale space scale invariant?

texts:
I'll use the Scale-invariant feature transform algorithm as an example here.  SIFT creates a scale space based on scaled gaussian filtering of an image, and then computes the difference of gaussians to detect potential interest points.  These points are defined as the local minima and maxima across the difference of gaussians.

It is claimed that this approach is scale invariant (among other puzzling invariances).  Why is this?  It's unclear to me why this is the case.


texts: noise remains after applying median filter

texts:
So I'm apply the median filter to an image, but at the output, there's blue semi-dots appearing. What are they?

This is the output

And the python code for anyone interested

from scipy.ndimage.filters import median_filter
from matplotlib.pylab import imread, imshow, subplot, title, show, get_cmap

img = imread('images/speckle.gif')
new = median_filter(img, 3)

subplot(121); imshow(img, cmap=get_cmap('gray'));
title('The original')
subplot(122); imshow(new, cmap=get_cmap('gray'));
title('Filtered with the median_filter')

show()


texts: Why Is The Maximum Likelihood Estimation (MLE) Method Taken as the Benchmark for Comparing with Other Methods?

texts:
In many research articles the performance of an estimation method is compared to that of the ML estimation performance. If the performance of the method does not achieve the ML estimation performance, then the method is 'suboptimal' or not good enough. I don't quite remember what is the meaning of ML estimation performance and sub-optimal or optimal. When is a method optimal and optimal with respect to what? I am following the book titled: "Statistical signal processing Vol 1 by S. Kay". As a beginner and self-learner it is quite hard to grasp the actual implication of the terms -- achieving ML estimation performance, optimal and sub-optimal.

Can somebody please explain an intuitive meaning to these terms? Thank you.


texts:
There have been many answered questions regarding time-invariance of a system, but none could
help me understand some cases I have had as an question:
$y(n) = 3x(n) + y(n-1)$
if the other answers in where to be trusted:
$x(n) = x(n+N)$ and $y(n) = y(n+N)$ should be equal:

$y(n+N) - y(n+N-1) = 3x(n+N)$

$y(n) - y(n-1) = 3x(n+N)$

which would mean it is time-invariant
but how about those 2 Equations if i do it the same why like i did before ill I still have the $0.9^n$
which is dependent on $n$ so it should be time-variant but how would i write this down the notation looks nice?
$y(n) = 0.8x(n-2) + 0.9^ny(n-3)$
$y(n) = 0.9a(n)x(n)$


texts: I have stumbled upon this term in many papers, though no one bothers to elaborate. Is this the case where I do not have prior knowledge of the array structure? or otherwise is it the number of microphones available? Maybe both?


texts: What is a good distance measure for matching SIFT descriptors depending on the distribution of their noise?

texts:
I have read some papers about distance measures like Euclidean, Manhattan or Chi-Square for matching gradient based image descriptors like those computed from the SIFT Algorithm (128-D vectors). Most of them state, that one or the other measure is more suitable for that task and that this decision depends on the noise distribution assumed (Like L2 is suitable if Gaussian noise is assumed or L1 for Laplacian noise).

So my question is: how do I determine the distribution of my descriptors and how does that affect the choice for a distance measure?

In one paper the authors stated that the distribution would be the distances between SIFT descriptors of correctly matching keypoints. But doesn't that mean this distribution is biased by the choice of the distance metric itself?!


texts: Why can't a causal digital filter have an infinitely sharp transition between the passband and the stopband?

texts:
In DSP book by Proakis and as well as in this pdf, it is mentioned that practical causal digital filters cannot have an infinitely sharp transition from Pass-band to Stop-band. Why is it so? Can you please provide a detailed explanation (with proof)?
Edit 1: In the book in was just mentioned that is a consequence of the Gibbs phenomenon, which results from the truncation of h(n) to achieve causality. I didn't understand how.


texts: What is the relationship between luminosity of Lab color space and intensity of HSI color space?

texts: I am wondering whether one can simply convert the intensity of a color in HSI space to the luminosity value in Lab space. I am referring to te CIELab color space as discussed here: https://en.wikipedia.org/wiki/Lab_color_space#CIELAB-CIEXYZ_conversions


texts:
I am performing convolution in MATLAB using the built in conv function. I am convolving an impulse response (length 500 samples)  with an 1800 sample long input signal. This input signal contains 450 NaNs at random indexes between sample 530 and 1285. The output signal of this convolution is length 2299, and now consists of 1255 NaNs.
I'm unsure why this happens, and how to avoid it?

I'm also not sure I've provided enough info, but if I can add anything else please let me know.


texts: How to design a good IIR smoothing filter?

texts:
I need to design a smooting filter for computing the envelope of a signal. My specs in Matlab are
fs      = 44100;                    % Sampling frequency of 44100 Hz
fpass   = 2;                        % Passband edge frequency of 2 Hz
Wp      = fpass / (fs/2);           % Calculate normalized passband frequency
fstop   = 5;                        % Stopband edge frequency of 5 Hz
Ws      = fstop / (fs/2);           % Calculate normalized stopband frequency
Rp_lin  = 0.0001;                   % Max. passband ripple of 0.0001
Rp      = 20 * log10(1+Rp_lin);     % Calculate stopband ripple in dB
Rs      = 10;                       % Min. stopband attenuation of 10dB

%%                      Design the elliptic lowpass filter

% Calculate the order of the elliptic lowpass filter

[N, W] = ellipord(Wp, Ws, Rp, Rs);

% Calculate the coefficients of the elliptic lowpass filter via Bilinear Transformation

[b_ellip_lp, a_ellip_lp] = ellip(N, Rp, Rs, W, 'low');

% Convert the filter coefficients to second-order sections

[sos_ellip_lp, g_ellip_lp] = tf2sos(b_ellip_lp, a_ellip_lp);

%% Distribute the gain equally over the sections
[N_sos, ~] = size(sos_ellip_lp);

sos_ellip_lp(:,1:3) = sos_ellip_lp(:,1:3) * (g_ellip_lp^(1/N_sos));

%% Normalize the coefficients

sos_ellip_lp = sos_ellip_lp / max(max(abs(sos_ellip_lp)));

%% Plot the 16-bit coefficients of the Elliptic low pass

% Quantize the coefficients for plotting

sos_ellip_lp_16bit = round(sos_ellip_lp * 2^(16-1));

figure(1)
freqz(sos_ellip_lp_16bit);
title('16-bit Elliptic low-pass filter');

I am designing an Elliptic filter with ellip. At first the gain is massive

After gain distribution over all biquads, it looks better.

After quantizing to 16 bit, I get a 10 dB gain!

I mean, I can just apply gain correction by trial and error. The filter is not adaptive and thus I can work with it. I would nevertheless like to not have these gain jumps. Filters for higher frequencies did not seem to have this issue.


texts: How to measure the distance between two circle-like-lines over 360 degree

texts:
I am looking for a slim and with that quickly programmable image processing strategy to measure the distance between two circles over 360°. It is important to notice that the circles are not perfect, see Fig. 1.

Fig. 1: Sketch of the circles

The application of this, so you can imagine the question a bit better, is: I have a cylindrical specimen with an radius r1(alpha) over the 360°, lets call the angle alpha. In the sketch it is the black circle.
Now I would like to measure the thikness "delta r" over alpha of a slim layer on this specimen. In the sketch the thin layer is painted in blue.

I draw the circles in kind of paint style on purpose as in the real image I would draw the line where I think both rings would be and would then start with calculating/image processing. I know that it would be possible to detect the edge with some kind of edge detection algo. (eg. Canny), but I am looking for a slim and fast way of getting to my result.

The perfect result could look like this: A function which is an approximation of the black line over alpha and then another function for the blue line.

Here a slim idea I had: Drawing two perfect circles over both cirles and substract them -> delta r, but the made error might be big.

So my question (1) is:
How would you solve this to get delta r over alpha? Thinking outside the box is welcome of course!

Question (2) is:
I know that the circles I drew are - from a mathematical point of view - not described at each alpha -> so I guess an approximation of the circle is somehow needed. Any ideas how to approximate?

Edit 1: Here you can see a real pictures. For now I don't have the optical system set up, so I just took a photo with my smartphone. Assume the length of the picture in x direction as 40 mm so you have a reference.

Fig. 2: Picture of the real specimen.


texts: Fast Fourier Transform for varying period (Order analysis)

texts:
I am trying to do a fast Fourier transformation on accelerometer data from a shaft rotating at varying speed in Python.

What I have done so far:

1: The original plot was in the time domain, and I therefore did a order analysis (resampled), and got the following plot:

This plot shows the angular rotation plotted against amplitude.

2: Now, an FFT was done with this code:

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

class FastFourierTransform:
    # Amplitudes is a row vector
    def __init__(self, amplitudes, t):
        self.s = amplitudes
        self.t = t

    # Plotting in the input domain before fft
    def plot_input(self):
        plt.ylabel("Amplitude")
        plt.xlabel("Shaft angle [Radians]")
        plt.plot(self.t, self.s)
        plt.margins(0)
        plt.show()

        '''
        The second half of this array of fft sequence have similar frequencies
        since the frequency is the absolute value of this value.
        '''
    def fft_transform(self):
        mean_amplitude = np.mean(self.s)
        self.s = self.s - mean_amplitude # Centering around 0
        fft = np.fft.fft(self.s)

        # We now have the fft for every timestep in out plot.

        # T is the sample frequency in the data set
        T = self.t[1] - self.t[0] # This is true when the period between each sample in the time waveform is equal
        N = self.s.size  # size of the amplitude vector
        f = np.linspace(0, 1 / T, N, )  # start, stop, number of. 1 / T = frequency is the bigges freq
        plt.ylabel("Amplitude")
        plt.xlabel("Frequency [Hz]")
        y = np.abs(fft)[:N // 2] * 1 /N

        # Cutting away half of the fft frequencies.

        sns.lineplot(f[:N // 2], y) # N // 2 is normalizing it
        plt.margins(0)
        plt.show()
        time = f[:N // 2]
        return fft, time

3. Result, with the normalized amplitudes plotted:

Questions:

Does this thought process look correct?
Is it correct to say that the final fft plot is in the frequency domain? From this link, http://zone.ni.com/reference/en-XX/help/372416L-01/svtconcepts/svcompfftorder/, it looks like the final plot domain should be in the order domain, but I am not sure since the fft was
done from the radian domain.

Vibration- Based Condition Monitoring of Wind Turbines by Tomasz Barszcz has this image

Thanks in advance.


texts: Graphing the response of an envelope follower

texts:
First off I want to apologize for not being very clear or technical in describing this question.  I haven't had any formal DSP schooling.
It will probably be easier to explain what I'm trying to accomplish.  In an audio plug-in I have a graph that displays the envelope used in a compressor.  To do this currently, I create an envelope detector with the same attack and release values as the compressor and then feed it samples of 1 until it reaches the full attack time, and then feed it 0s until it almost fully decays to 0.  I store the output at each sample point and use that to draw a plot.  What I would like to do instead is just figure out how to graph this so that I can make use of smooth bezier curves.
Here is what it looks like currently

The detector uses a discrete function to calculate the envelope based on the gain/release constants and the previous envelope value.
T detect(T input)
{
    T envIn = processInput(input);

    if(_envelope < envIn)
        _envelope = _envelope * _ga + (1 - _ga) * envIn;
    else
        _envelope = _envelope * _gr + (1 - _gr) * envIn;
    return _envelope;
}

Is there a way to transform this from a discrete equation into something that can easily be graphed?
Again I apologize for not knowing how to describe this in a mathmatical way, and I hope I'm not out of place for asking for help on this.


texts: Gaussian noise with different SNR levels are usually used in research works to simulate a realistic environment. How can researchers guarantee that Gaussian noise can simulate the reality of a System?


texts: Covariance matrix $P$ and $Q$ for complex valued data

texts:
The equations for the Kalman filtering are

  IM = H*X;
  IS = (R + H*P*H');
  K = P*H'/IS;
  X = X + K * (y-IM);
  P = P - K*IS*K';

The covariance matrix in my application is of dimension 3 by 3 size because the there are 3 state variables that need to be estimated.

For data in the real domain, I initialized the covariance matrices $Q$ (process covaraince) and $P$ (prediction covariace) in Matlab as

P= 5*(eye(3,3))
Q  = 0.001*(eye(3,3));

For the complex domain, I initialized $P$ and $Q$ as

    P_real = 5*(eye(d,d));
    P_imag = 5*(eye(d,d));
    Q_real = eye(d,d)* 10^(-2);
    Q_imag = eye(d,d)* 10^(-2);

    P  = P_real + sqrt(-1)*P_imag;

    Q  = Q_real + sqrt(-1)*Q_imag;

R is a scalar real valued number.

Problem  : The default operation in Matlab is complex conjugate represented by the operator (.) '. I am using a toolbox which applies the above expression for Kalman filtering and there is no term conj() present especially when evaluating the expression IS , K, P. However, the toolbox never uses any (.).' or (.).* operators as well. When in complex domain, would the equations change to

IM = H*X  % no conjugate

IS = (R + conj(H)*P*ctranspose(H))

K = P*ctranspose(H)/IS;

X = X + K *(y - IM)  % here * denotes the multiplication operation

P = P - conj(K)*IS*ctraspose(K)

where H, X, y, P are all complex valued.

Are these okay?


texts: How to Sample Smooth Random Walks of Instantaneous Frequency

texts:
I want to generate a random sequence for instantaneous frequency with as little aliasing as possible (i.e I expect to see as close to constant amplitude as possible, with only the frequency varying). What I have so far seems close, but in some cases results in a few frequencies with less-than-expected amplitudes:

def ms(x):
    return (np.abs(x)**2.0).mean()

def normalize(y, x=None):
    if x is not None:
        x = ms(x)
    else:
        x = 1.0
    return y * np.sqrt( x / ms(y) )

def brown(N, state=None):
    state = np.random.RandomState() if state is None else state
    uneven = N%2
    X = state.randn(N//2+1+uneven) + 1j * state.randn(N//2+1+uneven)
    S = (np.arange(len(X))+1)# Filter
    y = (np.fft.irfft(X/S)).real
    if uneven:
        y = y[:-1]
    return normalize(y)

inst_freq = np.cumsum(brown(100))
inst_freq -= inst_freq.min() - 1
inst_freq /= inst_freq.max()
rectified_inst_freq = 1.0 / (np.floor(1.0 / inst_freq))

plt.plot(inst_freq)
plt.plot(rectified_inst_freq, 'r')
plt.show()
plt.plot(np.cos(np.cumsum(np.pi * inst_freq)))
plt.show()
plt.plot(np.cos(np.cumsum(np.pi * rectified_inst_freq)))
plt.show()

The rectified_inst_freq variable is attempting to round the instantaneous frequencies to the a frequency with an integer period, this helps a lot with the aliasing, but not entirely. Results that show the strange high-frequency amplitude problem (the middle graph before rectification, and the bottom graph after):

The high frequency portion at the end of the rectified signal should have amplitude 1, but it is stuck at .75. I think this has to do with aliasing, but I am not sure how to fix it. In some scenarios (particularly ones that start at a high frequency) the amplitude is correct (see below)


texts: What is the basic concept of Fourier transform

texts:
I use discrete Fourier transform for digital image processing purposes , but I don't understand basic concept behind it. For example :

What information exists in frequency domain?
What is difference between spatial domain and frequency domain?


texts: Intuition behind the scaling property of Fourier Transforms

texts:
The Fourier transform of $f(ax)$ is $\frac{1}{|a|}F(\frac{u}{|a|})$. So the frequencies are scaled horizontally but the magnitudes are also scaled when the graph of $f$ is scaled horizontally.

On the other hand the fourier series of $f(\theta)$ on $[-\pi,\pi]$ is $\sum_{-\infty}^{\infty} c_n e^{in\theta}$ whereas the fourier series of $f(a\theta)$ on $[\frac{-\pi}{a},\frac{\pi}{a}]$ is $\sum_{-\infty}^{\infty} c_n e^{ina\theta}$ so that the frequencies are  scaled horizontally but the magnitudes are left unchanged at $c_n$.

Now these are easy to show mathematically, but I'm left feeling uneasy conceptually. Is there same intuition geometrical view to leave my intuition satisfied that explains why there is this lack of symmetry?


texts: Calculating frequency and damping ratio from transfer function given eigenvalues

texts:
I have the following standard transfer function for a damped linear oscillator:

$$G(s) = \dfrac{\omega_0^2}{s^2 + 2\zeta\omega_0s + \omega_0^2}$$

Now I have two eigen values at locations $-100 \pm 100i$.
I want to calculate the damping coefficient, $\zeta$, and natural frequency $\omega_0$.

I thought that the characteristic polynomial of matrix $A$ from the linear system in state space form is the denominator of $G(s)$. However, when solving $s^2 + 2\zeta\omega_0s + \omega_0^2=0$ for $s = -100 - 100i$, I get an equation with two unknowns right? What theoretical point am I missing here? how should I go about this?

Thanks in advance


texts: Why do object-boundaries look better on a badly rectified image?(openCV stereo correspondence)

texts:
I got 2 pictures which is almost parallel, and not positioned very far from each other.
I'm using OpenCV to try to create a disparity map (Stereo correspondence).
Because I'm trying to use it in a real world scenario, the use of chessboard calibration is a bit un-practical.
Because of that, I'm using stereoRectifyUncalibrated().

I tried to compare the results, using 2 different sets of corresponding points for the rectification:

Points manually selected(point & click)
Points generated from SURF and filtered with RANSAC

Input image1:

Input image2:

(Note that I do undistortion on the images before using them for rectification etc)

Rectified images with SURF and RANSAC:(1 and 2 in that order):

Rectified images using the manually selected points(which is more inaccurate!):

Now, the thing is, looking at the result we see that the surf-version is almost perfectly rectified.(The epipolar lines are quite well aligned).
While the manually selected point version is quite badly rectified...the epipolar lines are nowhere near aligned.

But when we look at the result of openCV's sgBM() using both our rectifications:

Manual point result:

SURF point result:

The disparity/depth shown is more accurate/correct with the SURF-point(well rectified version). No surprise there.
However, the actual detected object-pixels and object-boundaries are actually a lot better on the badly rectified verison.
For instance, you can see that the pen is actually a pen and has the shape of a pen, in the bad rectified disparity map, but not in the well-rectified map.
 Question is, why?
And how can I fix it?
(I tried fooling around with the sgBM() paramteres, making sure they are the same for both etc, but it does not have any effect. It is the different rectifications alone that makes the badly rectified image look good for some reason(with respect to object-boundaries)).


texts: Blind Estimation of Signal Parameters for Zero mean Normal distribution

texts:
Let $\mathbb{s}$ be zero mean with a normal distribution. Let one of the observed set of signals be given by $K\mathbb{s}$ and another by $L\mathbb{s^*}$, where $K$ and $L$ both are constants, and $\mathbb{s^*}$ is defined by complex conjugate of $\mathbb{s}$. Also the values of $KL$, $\mid K\mid$, $\mid L\mid$ are known, along with the known relation that $\mid K\mid^2 + \mid L\mid^2=1$.

Is there a way to estimate parameters $K$ and $L$? I tried solving non-linear set of equations, splitting real and imaginary terms of products $K\mathbb{s}$, $L\mathbb{s^*}$, however the system of equations seems to be underdetermined.

The set of equations by separating real and imaginary parts, I tried to solve were as follows:
$$
K_{re}s_{re}-K_{im}s_{im} = A\\
K_{re}s_{im}+K_{im}s_{re} = B\\
L_{re}s_{re}+L_{im}s_{im} = C\\
-L_{re}s_{im}+L_{im}s_{re} = D\\
K_{re}L_{re} - K_{im}L_{im} = E\\
K_{re}L_{im} + K_{im}L_{re} = F
$$
where $A$ & $B$ are the real & imaginary parts of $K\mathbb{s}$ respectively,  $C$ & $D$ are the real & imaginary parts of $K\mathbb{s}^*$ respetively and  $E$ & $F$ are the real & imaginary parts of $KL$ respectively.

After substitutions, I had following set of equations:
$$
s_{re}=\frac{CL_{re}+DL_{im}}{|L|^2}=\frac{AK_{re}+BK_{im}}{|K|^2}\\
s_{im}=\frac{CL_{im}-DL_{re}}{|L|^2}=\frac{BK_{re}-AK_{im}}{|K|^2}
$$
and got the relations between real & imaginary parts of $L$ and $K$ as follows:
$$
L_{re}=\frac{FK_{im}+EK_{re}}{|K|^2}\\
L_{im}=\frac{FK_{re}-EK_{im}}{|K|^2}
$$
However, I could not find a way forward to estimate $K$ and $L$, as further substitutions could not give me unique solutions for the parameters.
I hope someone has a better idea to attack the problem. Thanks in advance!


texts: Signal Processing using Fourier Transform

texts:
How can I derive the fourier transform of

                                   g1(t) = {A cos(2πfc t),−T/2 < t <T/2
                                            0,               elsewhere

and;

                                  g2(t) = cos(t), −∞ < t < ∞

my answer for g1 is;

∫ g(t)exp(-j*2*pi*f*t)dt = ∫ Acos(2*pi*fc*t)exp(-j*2*pi*f*t)dt

   =∫ 1/2A (exp(j*2*pi*fc*t)+exp(-j*2*pi*fc*t))exp(-j*2*pi*f*t)dt

   =∫1/2A(exp(j*2*pi*fc*t)exp(-j*2*pi*fc*t)+∫1/2A(exp        (-j*2*pi*fc*t)exp(-j*2*pi*fc*t))

   =1/2A∫exp(-j*2*pi(f-fc)t) + 1/2A ∫exp(-j*2*pi(f+fc)t)dt

   =1/2A(f-fc)+1/2A(f+fc)


texts: Frequency modulation: which algorithm for discrete signals (Java, Swift, etc.)

texts:
I'm writing an audio app in Swift - and one task consists in frequency modulation of a given signal (microphone input or sound file) - thus I'm iterating through all frames of a buffer like
let leftSourceData = readBuffer.floatChannelData?[0]
let rightSourceData = readBuffer.floatChannelData?[1]

let leftTargetData = renderBuffer.floatChannelData?[0]
let rightTargetData = renderBuffer.floatChannelData?[1]

for i in 0..<Int(readBuffer.frameLength) {

   // TODO: perform frequency modulation on source data, writing the results to target data
   leftTargetData[i] = ???
   rightTargetData[i] = ???

}

The equation for Frequency Modulation is clear - you'll find in the Internet s.th. like that:
 f(t) = A sin(2*pi*fc*t + I * sin(2*pi*fm*t))
 with fc = carrier frequency, fm = modulation frequency
 and I = modulation index.

The above-mentioned equation only covers a sinusoidal input signal - anyhow, I want to modulate an arbitrary signal - thus s.th. different than a pure sine wave.
For amplitude modulation, I've used the following algorithm (within the loop) which works well for Swift and iOS:
let val: Double = sin(Double(2 * modulationFrequency) * Double(index) * Double.pi / Double(renderBuffer.format.sampleRate))
leftTargetData?[i] = Float(val) * (leftSourceData?[i] ?? 0)
rightTargetData?[i] = Float(val) * (rightSourceData?[i] ?? 0)

How can I do the same for FM?


texts: Unexpected behaviour when using lfilter with initial filter delay multiple times

texts:
I want to filter a very long signal in smaller parts, therefor I am currently using scipys lfilter with an initial filter delay doing multiple iterations of:

(signal_part_filtered, filter_delay) = lfilter(b, a, signal_part, zi = filter_delay)

The first few iterarions are going as expected (original is blue; filtered is red):

But after a few iterations the inital filter delay keeps getting bigger:

I think the issue might be the limited precision of the float datatype, is there any way to solve this?


texts: Experimental determination of impulse response function

texts:
I am doing an experiment (I major in physics, not familiar with signal processing so there may be some naive mistakes, please be kind) which requires acting arbitrary pulses on a system.

I use an AWG to generate the pulse (current), bring it to the linear amplifier, and eventually, to a coil (than the current is converted to the magnetic field and act on the system). The pulse will be distorted by the coil and the amplifier, and U want to know the impulse response function (f), so as to modify my input pulse.

Is there a method which can determine the impulse response function experimentally?

I know a little about the method which says: generate a number of input pulses, use an oscilloscope to get the distorted pulse, and f is treated as a matrix which can be determined by multiply the "inverse of input pulse matrix" and the output pulse matrix. Can this really work? Is there a standard method to determine f (I suppose there is)?


texts: Need a 2 parameter ranking system for image compression

texts:
I'm working on a project in which specific bit planes of an image are extracted and compressed. This process is repeated for every combination of bit planes possible.For an image, I have calculated the compression ratio and the Structural Similarity Index(SSIM).

I need a method to allot ranks when compressed with different combinations of bit planes.

My progress:
I tried to give equal weight to the compression ratio and the SSIM value.

score = 0.5*ssim + 0.5*compression_ratio;

This equation gave me a value between 0 and 1. The problem with this equation is it gives an equal score for multiple combinations and few combinations are extreme opposites on the scale. For example, i got equal scores when I compressed an image using the {1,2,3,5,7,8} bit planes(High SSIM value) and only the 3rd bit plane(High compression ratio).

I need help in balancing the parameters and getting an appropriate score. The score doesn't need to give equal weights like I have given.

Thanks in advance


texts:
I am a little confused about definition of intensity scaling. It is something like histogram equalization or a method in which we just multiply all intensities by a scalar.
Thanks


texts:
So we're doing a project where we've decided on QAM-n modulation. I've read up a little bit on synchronization and found Barker sequences, which supposedly has nice autocorrelation properties which makes them good for synchronization.

So one option is to use BPSK modulation only in the synchronization phase and then move onto QAM for the data. But to me, it's very unintuitive that we couldn't construct even NICER synchronization sequences (wrt autocorrelation peakedness etc) if we utilize more of the signal space. Maybe something like simultaneous Barker BPSK for both I and Q?

Is there a nice explanation for why it's generally best to use BPSK for the synchronization?


texts:
If I have the histogram of an input image as Gaussian probability density function of the form:

$$P_r(r)=\dfrac{1}{\sqrt{2\pi}\sigma}e^{-\dfrac{(r-\mu)^2}{2\sigma^2}} $$

where: $\mu$ and $\sigma$ are the mean and standard deviation of the Gaussian PDF respectively.

The approach is to let $\mu$ and $\sigma$ be the measures of average grey level and contrast of a given image. What is the transformation function that can be used for histogram equalization?

Now, I found that I have to take small value for sigma and make some integration. So any help to do that please?


texts: Bandwidth required for orthogonal and antipodal signals

texts: What is the bandwidth required to transmit baseband modulated signals like an antipodal and orthogonal signals?


texts: What happens at the edges in fft2 filtering?

texts:
I am doing a fairly simple homework assignment in MATLAB that is to filter an image in the frequency domain. So my code looks something like this:

output = ifft2(fft2(image, x, y).*fft2(filter, x, y));

After testing some images, this code filters properly, but has weird edge behavior. Out of curiosity, what exactly happens at the edges using this method? (0, symmetric, replicate, circular)?


texts:
I have some noise measurements and 4 different models that can generate a similar noises using estimated parameters from the real measurements. To compare the estimated noises I calculated the PDF, the CCDF and the PSD.

To compare the results I used the Kull-Back Leibler divergence of the PDF, the MSE of the CCDF and the RMSE and the corr2 matlab function for the PSD.

I am not sure if I did the right thing.

Can you please tell me if what I did is correct and do you know any other way or function to make such comparison?

Thanks.


texts: Absolute maximum number of possible PN sequences for given length

texts:
I previously asked if it was possible to use encryption-grade randomness as a CDMA spreading sequence and learned that the codes must be orthogonal, not necessarily random.

I've tried different ways of generating codes and very few can coexist without mutual harm. Everything I have seen on the subject says that you can have up to n orthogonal possibilities for an n-bit spreading code. Is there any type of PN sequence that can produce more than n orthogonal CDMA codes?


texts: How to label sound events in noisy .wav files

texts:
I'm trying to find a way of automatically identifying birdsong in phone-recorded sound files (1min long). Currently the algorithm I'm using doesn't label all of the events I want it to. (it is designed for non noisy audio files).

I don't have much experience in DSP. I'm hoping that upon looking at my wave-forms someone may have an idea of how to pick out these 'birdsong events'!
I've shown one of the waveforms. The rest can be seen here: imgur.com

Any method would be appreciated!

I've thought about applying a noise reduction algorithm but I suspect this could just interfere with things.
Here is a link to the algorithm I'm using in case you're interested: Github/kylemcdonald

The labelled events are the cross hatched ones.

Thanks in advance for taking the time to look at this!


texts: Why is the size of results from FFT half the size of the input, while that is not the case in image processing?

texts:
When dealing with real data, Fourier Transform produces symmetric results and FFT algorithms discards half of it. An input with 1024 points will yield only 513 points in the output.

However, in image processing, when Fourier Transform is calculated, the size of the output image is same as that of the input image. Assume that the images are 2D. The (pixel) value of image function f(x, y) is still real. So, why do we still have an output image size same as input?


texts: Difference between single layer and dual layer beamforming in LTE

texts:
I'm trying to understand the difference between single and dual layer beamforming in LTE as described in TS 36.101 section B.4.1 and B.4.2 respectively. Also refer to TS 36.211 section 6.3.4.2.3 for selection of precoding matrices.

My questions are as follows,

What does the following linear operation(From TS 36.211 section 6.3.4.2.3) of subtracting twice the projection matrix from identity signify and how is it relevant to beamforming?
$$
\rm W_n=I-2u_n u_n^H/u_n^Hu_n
$$
What is the difference between single and dual layer beamforming based on their respective precoding transformations given below?

For single layer beamforming (From Section B.4.1 TS 36.101)
$$
\begin{bmatrix}y_{bf}(i)\\\tilde{y}_{bf}(i)\end{bmatrix}=\frac{1}{\sqrt 2}\left(W_1(i)y^{(7)}(i)+W_2(i)y^{(8)}(i)\right)
$$

For dual layer beamforming (From Section B.4.2 TS 36.101)
$$
\begin{bmatrix}y_{bf}(i)\\\tilde{y}_{bf}(i)\end{bmatrix}=W(i)\begin{bmatrix}y^{(7)}(i)\\y^{(8)}(i)\end{bmatrix}
$$


texts:
I wonder why several RMS meters I'm using for music production show RMS values above 0.

For example, when using an input signal with a constant DC of 1 (0 dbFS), I get these results, although Wikipedia says that in such cases, the RMS is equal to the amplitude of the direct current.

Samplitude Pro X (RMS: 2.9)

TT Dynamic Range Meter (RMS 3.0)

Whereas other meters show the values I expect:

Voxengo SPAN (RMS: 0.0)

Why is it that meters come to such different results?


texts: In OFDM, does the parameter time offset is important if there is only one mobile in the cell?

texts:
I am learning about OFDM, and in my reading, I have noticed that (mentioned everywhere) synchronization is a very critical aspect of OFDM. Here, I am writing in the context of time synchronization (time offset).
What I understood that time offset is important if there are multiple transmitting mobiles in a cell, and due to different delays in each path it is difficult to find the frame starting time at the receiver. Will this concept apply if there is only a single mobile transmitting in a cell?


texts: Proving that the IDTFT is the inverse of the DTFT?

texts:
The DTFT is given by:

$$X(e^{j\omega}) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}$$

The IDTFT is given by:

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

I have been able to show by substitution of the DTFT into the IDTFT that the transform and a subsequent inverse transform return $x[n]$:

$$\begin{align}
x[n]&=\frac{1}{2\pi}\int_{0}^{2\pi}X(e^{j\omega})e^{j\omega n}d\omega\\
&=\frac{1}{2\pi}\int_{0}^{2\pi} \left( \sum_{k=-\infty}^{\infty}x[k]e^{-j\omega k} \right)e^{j\omega n}d\omega\\
\end{align}$$

Swap the order of integration and summation:

$$x[n]=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\int_{0}^{2\pi}x[k]e^{j\omega (n-k)}d\omega$$

Argue that given $e^{j\omega (n-k)}$ is an odd function, it will only evaluate to anything other than 0 when $k=n$:

$$\begin{align}
x[n]&=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\int_{0}^{2\pi}x[k]e^{j\omega (n-k)} d\omega \ \delta[n-k]\\
&=\frac{1}{2\pi}\int_{0}^{2\pi}x[n]d\omega \\
&=\frac{2\pi}{2\pi}x[n]\\
\end{align}$$

However, I have been unable to show the dual case: that the inverse transform (IDTFT) substituted into the forward transform (DTFT) gives $X(e^{j\omega})$. How can we show this?


texts:
I was told that this question is too broad. Maybe I don't understand what this forum is about. My question is basically the same. I am looking for references to range Doppler algorithms so I can learn how they work. I said monostatic radar because it is probably the most basic radar system. If someone can help, fine. Else, I'll keep looking on my own. Still very appreciative.

I am looking for radar range-Doppler algorithms, especially for monostatic radars, used to create range-Doppler maps. I am looking into books and papers, and I checked the internet. There is much talk about RDA (Radar-Doppler algorithm) and RDM (Radar-Doppler map) but I cannot find any documents or code which are tutorial in nature.


texts: Force Linear Phase for a FIR Filter Synthesized Using Berchin's FDLS?

texts: As a follow-on to this post, how would you force linear phase for a FIR filter you synthesized using Berchin's FDLS?


texts:
Will someone be able to explain short time energy function for me in a simple way explaining all the variables in the formula?

The issue that I have is in the formula the summation varies m like,

m=n-(N-1) to n

from my research what I understood is that n is the center sample of the window, so if the center of the window is, lets say 80 and if the window size is 160. Then there is a situation where “m” get a negative value discontinuing a loop in a programming environment. Is my understanding incorrect or is there a reason? Thanks in advance.


texts: solving recurrence relation with Z-transform (LTI and causal)

texts:
I just wanted to doublecheck answers for my sanity's sake (exam next week)

problem statement

recurrence relation, solve it

$y[n+1]= 35 + y[n]*0.5$

according to my teacher it will be such that the input sequence must be causal, so that we have re-arranged form such as

$y[n+1]-  y[n]*0.5= 35u[n]$
,where $y[0]=20$, and  causal input

My understanding is that we should go as follows

get $X(z)$
find what is $Y(z)$
then we can maybe do inverse Z-transform to get the y[k] output sequence in explicit form sequence
typically the most difficult phase seems to be the inverse Z transform so it can be tricky, but we have a table, and our teacher taught us to either use

z-table
long division
partial fraction expansion and z-table

Firstly, Z-transform both sides

$Z(y[n+1])-0.5Z(y[n])=35*Z(u[n])$

according to the z-table we receive

$z*Y(z)-z*y[0]-0.5*Y(z)=35*\frac{z}{z-1}$

combine  terms and re-factor and emplace y[0]=20

$Y(z)(z-0.5)=35*\frac{z}{z-1}+z*20$

$Y(z)=35*\frac{z}{(z-1)(z-0.5)}+20*\frac{z}{z-0.5}$

We designate R(z) for the purposes of partial fraction expansion, because it is not found in the z-table

$R(z)=\frac{z}{(z-1)(z-0.5)}$

however the rightsided term is found directly in the z-table
$Z^{-1}(20*\frac{z}{z-0.5}) = 20* (0.5^k)$

start partial fraction expansions
$R(z)/z = \frac{1}{(z-1)(z-0.5)}=\frac{a}{(z-1)}+\frac{b}{(z-0.5)}$, then we find that, $a=2,b=-2$

$R(z)=2*\frac{z}{(z-1)}-2*\frac{z}{(z-0.5)}$

$Y(z)=35*(2*\frac{z}{z-1}-2*\frac{z}{z-0.5})+20*(\frac{z}{z-0.5})$

apply inverse Z transform

$y[k]=20*(0.5)^k+70*u[k]-70*0.5^k$

$y[k]=u[k]*(70-50*0.5^k)$

tabulating some results

k=-1; y[-1]=0
k=0; y[0]= 20
k=1; y1= 45
k=2; y[2]= 115/2


texts: Feedback Filtered-x LMS algorithm: question about theory

texts:
I was reading some of papers about the Active Noise Cancellation; in particular about the Filtered-x LMS algorithm, also known as FxLMS.

It seems that classic FxLMS (also called Feedforward FxLMS) algorithm is suitable to cancel broadband noise (white noise), while its Feedback version (the Feedback FxLMS, also abbreviated to FB-FxLMS) is suitable to cancel a narrowband noise (like sinusoid + white noise), but not a broadband noise.

I read that this is due to the fact that the Feedforward FxLMS has the primary noise provided as reference signal, while the FB-FxLMS has not the reference signal initially available, but it needs to estimate it (using the error signal).
And this estimation process is possible only if the primary noise has a predictable nature, so it can't be a white noise.

Is there someone (maybe who has already worked with the Feedback FxLMS algorithm) who can confirm this property of the Feedback FxLMS?

PS: As requested in the comments, here are the block diagrams of the two algorithms:

FxLMS block diagram:

Feedback FxLMS block diagram


texts:
I am currently attempting to demodulate an audio file using MATLAB. After taking the fft of the signal, I receive this:

After some more coding, I found the following values:

Fs = 44100;
Window size = 441000 (by taking length(signal))
Fc = 11416

My guess is that this signal is a SSB with a full carrier, and I have read that to recover the original signal from the IF SSB signal, the SSB must be frequency shifted down to its original range of baseband frequencies.

The question here is: How could I shift the frequency spectrum of this signal to center at 0 Hz?

I have been reading up on the fft size/bins and its relationship to the sampling frequency. Will this value be needed? I have also tried using fftshift(), but to no avail. The goal here is to shift the frequencies to be centered at 0 Hz, use a filter to get rid of some of the frequencies, then use an ssb demod function to hopefully get the correct sound from it.

Thanks in advance, I appreciate you taking the time to help another guy out.


texts:
This is the spectrum of the signal:

And this is the signal:

On the FFT spectrum you can notice the high amplitude peaks and I think these frequencies are the main ones in the signal:

I think that the "low" harmonic must be clean, but I have the problem that low amplitude components are placed in every frequency interval.  The main problem is that I haven't got any model of the signal for a hypothesis about noise, and as result, I must detect noise components.  After this step, I must design a filter or use another method of filtering for denoising the signal.  Which methods can I apply for (1) detecting noise and (2) filtering?


texts: Why is the Fourier transform of a Dirac comb a Dirac comb?

texts:
This doesn't make sense to me, because the Heisenberg inequality states that $\Delta t\Delta \omega$ ~ 1.

Therefore when you have something perfectly localized in time, you get something completely distributed in frequency. Hence the basic relationship $\mathfrak{F}\{\delta(t)\} = 1$ where $\mathfrak{F}$ is the Fourier transform operator.

But for the Dirac comb, applying the Fourier transform, you receive another Dirac comb. Intuitively, you should also get another line.

Why does this intuition fail?